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Preface 

This report contains the 1988 annual progress reports by the Postdoctoral 
Fellows of the Center for Turbulence Research. It is intended primarily as a 
contractor report to the National Aeronautics and Space Administration, Ames 
Research Center. A separate report entitled, “Studying turbulence using nu- 
merical simulation databases -II,” covering the 1988 Summer program activities 
was released earlier this year. 

The primary objective of the CTR is to stimulate and produce advances in 
physical understanding of turbulence, in turbulence modeling and simulation, 
and in turbulence control. The primary means by which CTR seeks to achieve 
these objectives is by bringing together key individuals in fields bearing on tur- 
bulence to address diverse problems in turbulence. Postdoctoral Research Fel- 
lows and students conduct research in collaboration with staff of NAS A- Ames 
Research Center and Stanford faculty members, using a substantial array of 
research facilities provided by both institutions. 

Four thrust areas have been established for research: 

1. Fundamental modeling of turbulence 

2. Turbulence structure and control 

3. Transition and turbulence in high-speed compressible flows 

4. Turbulent reacting flows 

These program areas have been used in selecting the Research Fellows and also 
as a guide for turbulence research by graduate students. The CTR roster for 1988 
is provided in the Appendix. All Fellows with tenure of more than two months at 
the Center provided a written report outlining their study and accomplishments 
which appear in the following pages. The reports are grouped in the general 
areas of modeling, theory and simulation, compressible and reacting flows, and 
experimental research. 


Parviz Moin 
William C. Reynolds 
John Kim 
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The length scale for sub-grid-scale 
parameterization with anisotropic resolution 

By D. K. LILLY 1 


Abstract 

Use of the Smagorinsky eddy-viscosity formulation and related schemes for 
sub-grid-scale parameterization of large eddy simulation models requires speci- 
fication of a single length scale, earlier related by Lilly to the scale of filtering 
and/or numerical resolution. An anisotropic integration of the Kolmogoroff en- 
strophy spectrum allows generalization of that relationship to anisotropic res- 
olution. It is found that the Deardorff assumption is reasonably accurate for 
small anisotropies and can be simply improved for larger values. 

1. Introduction 

Numerical integration of the time-dependent equations of fluid dynamics to 
simulate evolution of the largest scales of motion has been applied for about 30 
years, initially in the field of weather prediction. Early workers were concerned 
over what to do about the scales of motion too small to be resolved by the 
computer, but the problem was often confused with the numerical errors intro- 
duced by finite difference algorithms. The technique that Smagorinsky (1963) 
introduced grew out of an empirical variable eddy viscosity method due, I be- 
lieve to R. Richtmyer (no reference available), for smoothing simulated shock 
wave calculations. Smagorinsky and his associates recognized, however, that 
Richtmyer’s variable viscosity was also consistent with the notion of a universal 
equilibrium range of turbulence. I pursued this point and used the Kolmogoroff 
inertial sub-range hypothesis to quantitatively relate the length scale required 
in the Smagorinsky formulation to the resolved scale (Lilly, 1966, 1967). I as- 
sumed isotropic resolution, although numerical simulations then and now often 
are carried out with anisotropic resolution. The question of how to generalize 
my expression for this purpose has engaged some attention. Deardorff (1970) 
assumed that the length scale is proportional to the cube root of the product 
of the resolution scales in the three directions. Alternatives have also been pro- 
posed, most of them involving the effects of a nearby boundary. Piomelli, et al 
(1987) review the subject and conclude that no real consensus has yet emerged, 
though Deardorff ’s formulation is widely used. 

1 Permanent Address: CIMMS, University of Oklahoma 
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The purpose of this note is to aid in resolving this issue. I extend my earlier 
calculation, oriented now toward spectral and pseudo-spectral integration tech- 
niques, to allow resolution to vary with direction. I still assume, however, that 
the small scale limit of the simulation is in the Kolmogoroff inertial sub-range in 
all directions. This is a serious limitation, as the most common reason for apply- 
ing anisotropic resolution is an expectation of anisotropic and/or inhomogeneous 
turbulence, typically in regions close to a boundary. Thus the scale anisotropy 
problem is usually compounded with a boundary layer problem and with the 
likelihood that a classical inertial sub-range may not exist, at least in the scales 
resolved by the simulation. Piomelli, et al (1987) discuss and test several addi- 
tional parameterizations intended to improve the boundary layer resolution, but 
all are modifications of an isotropic formulation related to Deardorff’s. The re- 
sults to be presented here are offered to test and perhaps improve the Deardorff 
assumption. 


2. The problem and the solution 

Consider the problem of integrating the Navier-Stokes equation under condi- 
tions of incompressibility and constant density, so that 


dui /dt + (9(uiti/) / dxj -f dp/dxi = ud 2 Ui jdx 2 (l) 


and 


du{/dxi= 0 , ( 2 ) 

where p is pressure divided by density, and other terminology is conventional. A 
low-pass spatial filter is applied to these equations, in recognition of the limited 
resolution available to any real discretization technique. In the early develop- 
ment era this was typically a uniform average over a box surrounding a spatial 
grid point. Leonard (1974) introduced the use of Gaussian filters, a technique 
which is widely applied in engineering fluid dynamics simulations. In many 
current research simulations the effective filter is the spectral cut-off of a finite 
Fourier transform. In any case it is assumed to be a linear operator and is here 
designated by angle brackets <>. Because it is linear it may be exchanged with 
the differential operators in (l) and (2) and applied to the flow variables directly. 
This allows the linear terms of the equations to be integrated as if they were 
unfiltered. The filtered momentum flux product, < Ujtiy >, cannot, however, be 
resolved directly from the filtered velocities, and must be subject to a creative 
parameterization . 

I assume the decomposition of the momentum flux product applied by Sma- 
gorinsky, Lilly, and Deardorff, i.e. 
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< U,Uy >—< Ui >< U y > —Tij + 

Tij = — < U,Uy > + < Ui >< Uj > —SijE 

E = (<u 2 > - < Ui > 2 )/2 (3) 

with Tij designated as the subgrid scale (SGS) stress, E the SGS kinetic energy, 
and 6ij is the Kronecker delta. The filtered equations of motion and continuity 
are now written as 

d < Ui > /dt + 3(< Ui >< Uj >)/dxj + dn/dxi 

= drij/dxj + ud 2 < Ui > /dx 2 (4) 

d < Ui > /dxi — 0, (5) 

where i r =< p > +E. 

The Smagorinsky parameterization for the SGS stress is of the eddy viscosity 
type, given by 

nj = K(Sij + Sji), ( 6 ) 

where S^y — d < Ui > /dxj, and 

K = A 2 5, (7) 

where S 2 = S»y(S,y + Sji). The length scale A was assumed by Smagorinsky to 
be proportional to the grid resolution. Some investigators have replaced S by 
lo, the square root of the enstrophy, i.e. to 2 = (3u, / <3xy — 3uy/3x t ) 2 / 2. Not 
much difference is found in practice, corresponding to the fact that the volume 
integrals of S 2 and w 2 are identical and they have the same Fourier spectral 
amplitudes. 

I evaluated the length scale A by assuming the validity of inertial cascade the- 
ory and neglecting intermittency effects. If the resolved kinetic energy equation 
is derived by multiplying (4) by < u,- >, the contribution from the eddy stress 
term is 

< Ui > drij/dxj = d(< Ui > Tij)/dxj — KS 2 

= d(< Ui > Tij)/dxj - A 2 S 3 (8) 

The first term on the rhs is variable in sign and vanishes or is small in the 
volume average, while the second term is always negative. It is assumed that 
the loss of energy from the resolved scales, A 2 S 3 , is the same as the viscous 
dissipation of energy, e. These are equated after evaluating 5 2 from spectral 
integration, i. e. 

S 2 = [ m k 2 E(k)dk 


( 9 ) 
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where E(k) is the spectral amplitude of kinetic energy, defined as an integral 
along the surfaces of spheres in spectral space. The wavenumber k m must cor- 
respond to the resolution limit, which is the cut-off wavenumber for spectral 
discretization. It is assumed that at and near k m the energy spectrum is given 
by the inertial range form, i.e. 

E(k) = ae 2 / 3 fc- 6 / 3 , (10) 

where a is the Kolmogoroff constant, about 1.5, and e is the viscous dissipation 
of kinetic energy. Substitution of (10) into (9) leads to the result of integration 

S 2 = (3/4 )ae 2 / 3 fc m 4 ' 3 (ll) 

The rate of energy removal from the resolved flow is evaluated by substituting 

(ll) into the last term of (8) so that 

A 2 S 3 = (3a/4) 3 / 2 A 2 ek m 2 (12) 

Equating this expression to dissipation allows evaluation of the length scale as 

A = S~ 3 /V/2 = ( 4 /3a) 3 / < ‘fe m - 1 (13) 

One may now set k m = constant/ A, where A is the grid interval, assumed 
to be isotropic. The minimum value for the constant, assuming the spectral 
wavelengths are bounded by the Nyquist limit, is two. When pseudo-spectral 
integration techniques are applied, in order to avoid aliasing the constant must 
be no less than three, which increases A by 50%. 

I now wish to drop the assumption of resolution isotropy. For simplicity 
I assume, however, axisymmetry, i.e. that the limits of resolution in two of 
the three dimensions are the same. The typical situation involves increased 
resolution close to a boundary, either through a variable grid length or use 
of Chebyshev polynomials. Within the axisymmetric framework one also may 
consider the effects of decreased resolution in one direction. In the conclusion 
section I argue that the results for different resolution in all three directions can 
be determined adequately from the axisymmetric case. 

Results are obtained by writing the enstrophy spectrum function in three- 
dimensional spectral space and integrating it over prolate or oblate spheroids, 
that is figures obtained by rotating an ellipse about its long or short axis. This 
may not correspond accurately to the resolution limits in a model using Fourier 
modes in two dimensions and Chebyshev polynomials or finite differences in the 
third, but it is more readily compared with the isotropic case. Similar results 
can be obtained by assuming a cylindrical volume in wavenumber space, but 
they differ only slightly and not in any qualitatively important manner. 
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If the three-dimensional spectrum function is isotropic, as assumed, it is 
E(k) / 4nk 2 , and the 3-d enstrophy or S 2 spectrum is E(k)/4ir. The wavenumber 
spheroid is assumed to be bounded by the equation for an ellipse, i.e. 

{kh/khm) 2 + ( k x /k zm ) 2 = 1 

where kh is the wavenumber in the axisymmetric direction, notationally assumed 
to be horizontal, and k x is that in the vertical direction, with kh m and k xm the 
horizontal and vertical spectral limits. In place of (9), the evaluation of S 2 is 
accomplished by integration over the spheroid as follows: 

S 2 = / / E{k)k h dk h dk x (14) 

Jo Jo 

The limit on the first integral is khi = ^ m (l — k x 2 /k xm 2 ) 1 / 2 - Note that if 
E{k) were replaced by 4n in (14), the integral would give the volume of the 
wavenumber spheroid, 47rfcj, m 2 fc zm . With the energy spectrum given by (10), 
with k 2 = kh 2 + k 2 , integration of (14) yields 

S 2 = (3/4 )ae 2/5 k hm 

where r = k hm /k xm and y = r 10 / 9 

The particular form for S is chosen so that if it is substituted into the first 
equality of (13), with y = 1 and the wavenumbers assumed to be inversely pro- 
portional to the grid spacings Ax, Ay, A z, one obtains the expression assumed 
by Deardorff, i.e. A ~ (AxAyAz) 1 / 3 . More generally, upon substituting S into 
(13) one obtains 

X = A/>y -3 / 4 , Xd = (4/3a) 3 ^ 4 fcfc m 2 ^ 3 k xm ~ 1 ^, (16) 

with Xjo equivalent to Deardorff ’s expression. 

The integral for y in (15) can be reduced to the forms 

y = r */ 9 {l-r 2 )- 1 ' 2 r^cosx^dx.tan*™ = (1 - r 2 ) 1 / a /r,for r < 1 (17a) 
Jo 

y = r 4 / 9 (r 2 - l)" 1 / 2 I*™ (cos x)~ 2 t z dx, sin x m = (r 2 - l) x / 2 /r,for r > 1 (176) 

Jo 

For the extreme cases r = 0, oo, the integration limits x m become 7r/2 in both 
expressions, and the integrals are evaluated in terms of tabulated Gamma func- 
tions. Thus for small r, S 2 ~ khm*/ 5 and A ~ khm 1 , while for large r, 
S 2 ~ kh m 1 / 3 *. 

zm and A ~ khm 1 '*k xm 3 ^ 4 - A somewhat more revealing in- 
terpretation can be obtained by differentiating A with respect to r. This yields 
the following: 


8/9 fc^ 4/S y(r) (15) 

[ [r 2 + (1 - r 2 )x 2 ] -6/6 dx. 
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— 1/3 for r C 1 (18a) 

r\~ 1 d\/dr = —(3 r/4y)dy/dr = 0 for r = 1 (186) 

5/12 for r > 1 (18c) 

d 2 \/d(lnr) 2 = (4/27)A r> at r = 1. (18 d) 

These show that the Deardorff expression is valid in the vicinity of r = 1, but 
that with increasing anisotropy in either direction A becomes larger than A d- 
Fig. 1 is a plot of A/ A £>. The data for the curve labelled “Length scale” were 
calculated to two-three digit accuracy by using the (I8d) for .05 < r < 20 
and matching that with an expansion of the expressions in (l7a,b) around their 
limiting values. The straight lines are plots of the asymptotic forms, that is 
y~ 3 / 4 for the y's given by (17a) for r C 1 and (17b) for r> 1. 

An accurate approximation for .02 < r < 50 is 

A/A 0 = r( 2 / 27 ) 1/, +r-( 2 / 27 ) 1/J -l (19) 

This is obtained as a solution to (l8d), but with the rhs replaced by (2/27) (A -f- 
A.d). This relation is also plotted in Fig. 1, and is imperceptibly different from 
the “Length scale” except at the most extreme anisotropies. 

3. Conclusions 

The results of the above calculation indicate that the Deardorff formulation 
is accurate to within 20% for .2 < r < 5. For anisotropies greater than that the 
length scale should be increased, with (19) sufficiently accurate for anisotropies 
less than a factor of 50 either way. Such large anisotropies are unlikely to produce 
accurate simulations in any case, and probably need to be enhanced by more 
sophisticated parameterizations than the Smagorinsky formulation. 

Upon consideration of the nature of these solutions, it seems evident that 
they can be extended to the case where the resolution is different in all three 
directions. Since equation (18) or (19) is symmetric in Inr, it doesn’t seem to 
matter much whether two wavenumber limits are larger and one smaller, or one 
larger and two smaller. Apparently the results are only sensitive to the largest 
ratio of the length scales, so presumably the existence of a third dimension with 
an intermediate resolution would have little effect. 
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FIGURE 1. The turbulence length scale, normalized by DeardorfTs (1970) 
assumption, and determined by solution of (13) and (15), for various ratios of 
the limiting wavenumbers, k^m/kzm- The straight lines area symptotic forms 
from (l7a,b) and the curve marked with + signs is an approximation given by 
(19). 
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Turbulence modeling 

By T. -H. SHIH 

1. Motivation and Objective 
Motivation 

In recent years codes that use the Navier-Stokes equations to compute aerody- 
namic flows have evolved from computing two-dimensional flows around simple 
airfoils to computing flows around full scale aircraft configurations. Most flows 
of engineering interest are turbulent and turbulence models are needed for their 
prediction. Yet, we know that present turbulence models are adequate only 
for simple flows and do poorly in complicated flows such as three-dimensional 
separation, large-scale unsteadiness, etc. The same progress that allowed the 
development of these aerodynamic codes, namely the introduction of supercom- 
puters, has allowed us to compute directly turbulent flows, albeit only for simple 
flows at moderate Reynolds numbers. These direct turbulence simulations pro- 
vide us with detailed data that experimentalists have not been able to measure. 
This work is motivated by the fact that data exists for developing better turbu- 
lence models and by the need for better models to compute flows of engineering 
interest. 

Objective 

The objective of this work is then to develop turbulence models for engineering 
applications. The model categories that show promise for immediate use are on 
the two-equation level and the Reynolds-stress level. We will make use of existing 
methodologies to develop models. The models will be tested using data from 
direct simulations, experiments and analysis. Specifically, our objectives are as 
follows: 

1. Examine the Reynolds stress budgets using direct simulation flow fields 
(Mansour et al 1988, Moin et al 1989). 

2. Use Rapid Distortion Theory to analytically study the effects of mean 
deformation on turbulence. In particular, examine the development of the rapid 
pressure-strain under rapid distortions. 

3. Compare existing models with data and theory. Develop models where 
needed using appropriate expansions and constraints. Test these new models 
using results from direct simulation, experiment and theory. 

4. Use the method of moment generating function to extend second order 
closures to higher order closures. We know that there exist a close connection 
between high-order moments and coherent structures. The moment generating 
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function approach should allow us to include the effect of coherent structures on 
the models. 

5. Extend turbulence modeling to compressible flows. 


2. Work Accomplished 

1. Numerical simulation of a three-dimensional boundary layer — AIAA pa- 
per no. 89-0373 (work supported in part by AFOSR) 

The objective of this simulation is to study the mechanics of three-dimensional 
boundary layers and develop improved models for their predictions. Three- 
dimensional effects were achieved by direct simulation of a fully developed tur- 
bulent channel flow subjected to transverse pressure gradient. To obtain a good 
statistical sample during the transient period, 14 computer runs were ensem- 
ble averaged. Each run started from a different realization of the channel flow 
(far apart in time). The simulation shows that, in agreement with experimental 
observations, the Reynolds stresses are reduced and that near the wall a lag 
develops between the stress and the strain rate. The reduction in the stress is 
due largely to a drop in the production rate and an increase in the dissipation 
rate. In the coming year, we will study the performance of existing second order 
closures in predicting this drop and introduce improvements where needed. 

2. k-e modeling 

k-e turbulence models are often used in computing engineering flows with 
moderate success. In general, the flow field is qualitatively well predicted, but 
quantitative agreement often falls short. In the case of Reynolds stress modeling 
the question of a length scale for an eddy viscosity does not arise, and the energy 
dissipation rate, e, is one of the unknowns of the problem. In this work an e 
equation has been developed following the methodology of Lumley (1978). We 
used it in connection with both two-equation modeling and Reynolds stress 
modeling. 

A two-equation model for use in low-Reynolds number flows has been tested 
against the channel data of Kim, Moin and Moser (1987) . The modeled equations 
are given by, 


c ,t 





k,k + vrSijSij — e 

J.fc 



,* 


+ Ci-urSijSij 
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o Jfc = 1 
— 1.3 
Ci= 1.48 
C 2 = 1.8 
e = € — (}> 
tj* = kjkj/ (2fc) 
i*r = Cpfpk 2 /*. 

C„= 0.09 

/„ = 1 - exp(-0.0115— ) 



We find that the model adequately predicts the mean velocity profile (see figure 
l) and the turbulent kinetic energy (see figure 2). 

3. Second order modeling of near-wall low-Revnolds number turbulence 

A set of second order closure models for low-Reynolds number turbulence 
has been developed for the simulation of wall bounded flows without using wall 
functions. The wall effect is built in the pressure-strain correlation term of the 
Reynolds stress equation and in the modeled terms of the dissipation rate equa- 
tion. We find that realizability is particularly important for modeling the near 
wall turbulence. The proposed models are particularly suitable for surfaces of 
arbitrary topology since they do not use the wall distance as a parameter. The 
models are tested by computing the fully developed channel flow. The full set of 
equations are used to compute the mean velocity, all the Reynolds stresses and 
the dissipation rate of the turbulent kinetic energy. We find reasonable agree- 
ment between the prediction of the mean profile and the data (see Figures 3) 
which is an indication that the shear stress is well predicted. However, the nor- 
mal stresses are not as well predicted. In particular the streamwise component 
is underpredicted (see figure 4) while the transverse component is overpredicted. 
The cause of this shortcoming is still being investigated. 

4. Second order modeling of a passive scalar in turbulent shear flows — AIAA 
paper no. 89-0607 

A model equation for the scalar dissipation rate was proposed using the ansatz 
that the ratio of mechanical time scale to scalar time scale has an equilibrium 
value. In addition a model for the pressure related terms in the scalar flux equa- 
tion was constructed based on consideration of realizability. The models were 
tested by comparison with experimental data for heated plane and axisymmetric 
jets. 

5. Rapid Distortion Theory and 2-D 2-C turbulence modeling 

Rapid distortion theory was used to analyze the development of the Reynolds 
stress in a 2-D 2-C homogeneous turbulence under mean irrotational strain and 
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mean shear. Here, 2-C refers to two-component turbulence. We found that the 
development of the Reynolds stress under 2-D 2C conditions is very different 
from its development under 3-D 3-C conditions. The findings can be summarized 
as follows: 

1. The mean shear or strain have no effects on the isotropy of the Reynolds 
stress if the initial field is isotropic. 

2. Initially anisotropic fields become isotropic under the influence of mean 
shear or strain. 

These findings are opposite to the finding in 3D 3C turbulence where mean 
shear or strain will drive isotropic turbulence away from isotropy. In 2D 2C 
turbulence the pressure strain drives the Reynolds stresses back to isotropy. 

A general 2D 2C model for the pressure-strain term, Tii, was constructed 
which contains only one undetermined coefficient C, 

2^ = Si, / "H ^CbjkbpjSpic + 2fl jkbik + 2n ik b k ] 

where, 

Si,- = \ (U,.i + Vi.,), n,,- = i {V,J - Vi.,) 

, u7u7 1 _ , 1 

bij -~n- 3 %. k ~ 2 U ’ U ’ 

Rapid distortion theory will be used to determine this coefficient. We postulate 
that the more general 3D 3C model should reduce to this model in the limit of 
2D 2C turbulence. 

3. Future Plans 

1. Examine the performance of existing second order closures in predicting 
three-dimensional turbulent boundary layers. Introduce improvements where 
needed. 

2. Use Rapid Distortion Theory to study the effects of mean deformation on 
turbulence. Extend the analysis developed for mean strain and shear to study 
the effect of rapid rotation on the turbulence. 

3. Use the method of moment generating function to extend second order 
closures to higher order closures. 

4. Use non-weighted ensemble averaging method (Shih et al, 1987) to extend 
the second order models to compressible flows. This averaging technique (as op- 
posed to Favre averaging) will allow us to apply all the incompressible modeling 
methodologies to modeling compressible flows. 
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Figure 1. Mean velocity profile in a channel as predicted by a two-equation 
model. 
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On curve and surface 
stretching in turbulent flow 

By N. ETEMADI 1 

Cocke (1969) proved that in incompressible, isotropic turbulence the average 
material line (material surface) elements increase in comparison with their initial 
values. We obtain rigorously, among other things, good estimates of how much 
they increase in terms of the eigenvalues of the Green deformation tensor. 

Introduction 

In the following note we will study deformation of material curves and surfaces 
convected in a turbulent flow. We will do this by looking at the notion of material 
line and surface elements as a means of generating hypotheses for the so called 
flow path ( Lagrangian ) of the motion. Then we will use these hypotheses to 
obtain upper and lower bounds for the evolution of the ensemble average of the 
arc length (surface area) of an arbitrary curve (surface) in time. 

For the definition of material line and surface elements and their historical 
background see Monin and Yaglom (1975). Cocke (1969) is the first who gen- 
erated these mathematical assumptions, see section 1.1, “implicitly” and gave 
a convincing proof of them for isotropic turbulence. Orszag (1970, 1977) takes 
these assumptions for granted, and by a variant of Cocke’s arguments obtains 
somewhat weaker results than Cocke’s, see Remark 1.2. Our work complements 
the work of Cocke. Namely, we will bring out these assumptions in section 1, 
and we will show, first, that once one accepts these assumptions then Cocke’s 
(1969) results can be improved to obtain “tight” upper and lower bounds for the 
ensemble average of material lines and surfaces. Then, in the remaining part of 
section 1 and in section 2 we carry on the results to arbitrary curves and surfaces, 
and in turn we also obtain upper bounds for moments of the dispersion between 
two points moving in the flow at any given time in terms of their separation at 
initial time. 

Throughout our work we will use x(a, t) as the Lagrangian representation of 
the flow, i.e. the trajectory followed by the particle which is at a at initial time 
to. We will assume that x is smooth enough in a space-time region and is, for 
fixed t, an invertible mapping so that our manipulations are legitimate. We 
will also use |t>| as the magnitude of an arbitrary vector t;, whose components, 
without danger of confusion, will be denoted by Vi,t/2, and t> 3 . 

1 University of Illinois at Chicago 

\ 
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1. Curve Stretching 

The motion of an infinitesimal material line 61 is governed by the equation 


d6l d6l ^ c , 

— = («-V)u or — = 


(i) 


with 61 (t 0 ) = S°l. Where ^ is the material derivative following the motion, see 
e.g. Monin and Yaglom (1975) section 24.5. It is easy to check that the above 
equation has the solution, 

3 

61 = ^2x <k 6°l k ;x(o,t 0 ) = a , (xij(a,t 0 )) = (fty) = I, (2) 

k= 1 


where I is the identity matrix. 

Let p(a, t) be the Lagrangian density. It is easy to show that 


det(nj(a, t )) 


p(o» to) £o 
p(o,t) ~ p t ’ 


( 3 ) 


where the left hand side is the determinant of the Jacobian matrix (or the 
Jacobian) of the transformation y = x(a, t),t > to, see Batchelor (1977) p.79. 
In particular, the Jacobian is one when the flow is incompressible. 

Next consider 


\6l\ 2 = 6 0 l T (x i<j ) T (x iJ )6°l = 6 0 l T W6°l. (4) 

Clearly W is a symmetric non-negative definite matrix. Since its determinant is 
(po/pt) 2 , consequently its eigenvalues, say ivi,W 2 ,vJs, are strictly positive and 
we have 

2 

det(W) = det(Xi,j) 2 = iv x w 2 w 3 = (5) 

Pt 

Let A = [a,]) be the unitary (rotation) matrix corresponding to diagonalization 
of W. Thus |i45°/| = |5°/|. We can rewrite, 

\6l\ 2 ( A6°l)f (. A6°l)l ( A6°l)% 

\6°l\ 2 “ Wl \A6°l \ 2 + W2 \A6°l \ 2 + W * \A6°l\ 2 ( 6 ) 

= (sin 2 6 cos 2 ip) w i + (sin 2 Osin 2 ip) + (cos 2 0)ws, 

4 fOl 

where 0, ip are the usual spherical coordinates of the unit vector 

Now we are in a position to take ensemble averages of both sides of this equa- 
tion. To do so, we need to make some assumptions about the joint probability 
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distribution of the variables involved; namely Wi’s, 0 and t p. This is where the 
physics of the problem come in. The following three assumptions have been 
extracted from the work of Cocke (1969) who proved them to be true for incom- 
pressible, isotropic turbulence. 

Assumptions; 

(1) Wi’s are independent of 0 and ip, 

(2) 0 and ip are uniformly distributed over the unit sphere, 

(3) Wi’s are identically distributed. 

Assumption one implies that the random orientation of the element 61 is un- 
correlated with its random deformations along the principal axis. Assumption 
two simply says that 61 is equally likely to be oriented in any direction. Finally, 
assumption three means that the random deformations of 61 along its principal 
axis have the same probability law. 

Proposition 1.1. Under Assumptions (l) and (2) we have, 

^(v^ + V^+V^) < (|^||) < ^<vW+ v/^+V^). (0 

!<(«>! + + «*)*} < (j^jj) < <(u>l + + U/s)*>. («) 


Proof. Since y/x is concave and the coefficients of Wi's add up to one, from 
(6) and (19) we obtain, 

l«l 


i«°/i 

Also it is clear from (6) that, 

i«i 


> (sin^dcos^xJj)y/uy[ + {sin^0sin^rp)y/w2 + [cos 2 6) y /ws. 


i«°'i 


< \sinOcosip\y/wl -I- \sin0sinxp\y/w2 + \cos0\y/w3. 


( 7 ) 


( 8 ) 


Next, use the joint density function of 0 and ip, i.e. j^sinO 0 < 0 < it, 0 < 
ip < 27T, to conclude that the average of the coefficients of tuj’s and their square 
roots, in (6), are | and respectively, and this will in turn give us (*). Since 
the coefficient of to, ’s in (6) are bounded by one we are only left to show the left 
hand side of (tt). To see this note that again by the fact that yfx is concave, 
from (6) and (19) we have, 

| 5 /| , . j.,(sin 2 0cos 2 tp)w 1 + (sin 2 Osin 2 ip) w-2 + (cos 2 0)w3, J 

iToTT = ( w i + w 2 + u>s) 2 p - — ; 

^ , . i . Ism0cos^i|iyi + \sin 0 sinip\w 2 + Icos^Ium. 

> [wi +W2 + t«3)>[ J 1 1 '■ ! ]. 


1* 


lt>l + lt>2 + U>3 


( 9 ) 
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Now the average of the coefficient of Wi's in the above line is | and we are 
through by Assumption (1)./// 

Proposition 1.2. Under Assumptions (1) and (2) and/or (l) and (3) we have, 


j( w i + u>f + w$) < ((i^jj) P } < + wj) fo 0 < p < 2, (t) 


= ^(u>i +™2 + W 3 ) for p = 2, 


(*) 


< ( (M r) < { ^1±A±A } to , < „ < OC, (m) 


^(log{wixV2W 3 )) < 


(iv) 


Proof. In light of what we have said in the proof of the above proposition, 
(u) is immediate under Assumptions (1) and (2). Now we can conclude («) by 
utilizing Assumptions (1) and (3) in (6) to obtain, 

( |j^j| 2 ) = (sin 2 0 cos 2 ip) (tv i) + (sin 2 9 sin 2 tp) (tv 2 ) + (cos 2 0)(w 3 ). 

= (sin 2 0co8 2 tp + sin 2 dsin 2 tp + cos 2 0)(tv i) = (u>i) (10) 

= + tt>2 + W 3 ). 

The left hand side of (t) and the right hand side of (m) follow from an inequal- 
ity like the one given in (7). y/x should be replaced by , and the direction 
of the inequality should be reversed due to convexity when p > 2. The right 
hand side of (t) is trivially true and the left hand side of (m) is a consequence 
of < ((|^r|j) p )» P > 1, see (18), and the right hand side of (*) for p = 1. 

Finally to obtain (iv) use the concave function log(x) rather than yjx in (7) 
under Assumption (l) and (2), and an argument similar to the one given at (10) 
under Assumption (1) and (3 )./// 

Remark 1.1. Note that the left hand side inequalities in the above propositions 
are strict unless W! = W 2 = W 3 . For all the concave functions involved are strictly 
increasing. Clearly this happens only when W is the identity matrix, meaning 
pure rotation. Furthermore, this has to be the case with probability one in order 
to have equality in the above propositions, which is certainly not an interesting 
case. 
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) > 1. 


for any p > 0. 

Proof. Cocke (1969) gives an a priori proof to the effect that the Assumption 
(1) to (3) are true for isotropic incompressible flow. This invokes the above 
proposition and together with Remark 1.1 and the fact that detfaj) = 1, see 
(3), we obtain the first inequality. The second one follows from (18), and the 
one we have just established as follows; 

<1^) = - exp W° 9 ]iwjF>i = ex p{p( l °9 > i- (ii) 

We could have also use d the arithmetic geometric mean inequality, (w[ + + 

u>5)/3 > {/(u>ittf2U>3) r , to achieve the same end./// 

Remark 1.2. The proof of the Assumptions (l) to (3) is the thrust of the work 
in Cocke (1969) in which he has also shown Proposition 1.2 (ii) and (iv), the 
above remark and corollary by using the same argument. Orszag (1970) takes 
these assumptions for granted, follows Cocke’s argument and obtains a weaker 
result, ( ) > 1, see the above corollary for p = 2. Note that this result does 
not imply that the average material line stretches. Now if we agree on all three 
assumptions, then it is trivial to show that (W) = 7(^»y); where 7 is the right 
hand side of (ii) in Proposition 1.2. In this connection see also Orszag (1977), 
p.240-241. 

Next we extend the above results to an arc length following the flow. The 
statement of the inequalities needed to carry this on can be found in the ap- 
pendix. Let C(s;to) : [a, 6] — ► R 3 be a parametric representation of a non- 
random curve at time to. Then C(s;f) = x(C(s;t 0 ),t) is the corresponding 
random curve , following the flow at time t. Let C to and C t , t > to be their arc 
lengths, respectively. For homogenous turbulence define, 


o(() = ( vW+^£+v^ ) 5 m - T 

< 'lit) = \f( 


W1+W2 + U >3 


) 


u>i + u;2 + W3 


) 


where the first and second inequalities are the consequence of (19) and (18) 
respectively. 

Theorem 1.1. For an isotropic, incompressible flow, 

Ct 0 < exp{(log( C t ))} < (C t ) < i(t)C to , 


(0 
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C t0 < a(t)C to < (C t > < ^a(«)C t0 , (it) 

^J-0(t)C to < (C t > < V3/?(t)C t0 . (Hi) 

Proof. The above corollary implies that for the non-random vector 6 °l = 

C"(s;« 0 ), 


(Ml 2>,*C£(s;t 0 )|)> = (Ml«l)> > log(\ 6 °l\) = M|C'(s;t 0 )|). (12) 

fc = l 


Consequently, 


exp{(/og(C t ))> = exp{(log[ j <*«]» 

(by (20), Si = [a, 6]) > J exp{(log (\ 9 - —^’ -|))} ds 

rb $ 

= exp{(log(\Y^ x ,k(C(s;t 0 ),t)C' k (s;t 0 )\)}} ds (13) 

Ja fc= l 

(by (12)) > [ exp{/op(|C'(s;t 0 )|)}ds 

J a 

= [ b \C'(s]t 0 )\ ds = C t0 . 

J a 


Next, inequality in (t) is an immediate consequence of (18) with <j> = exp(x). 
The following one is true by virtue of (18) with <f> = y/x, and Proposition 1.2 
(ii). For the rest of the inequalities all we need is Proposition 1.1./// 

Remark 1.3. Since in an inviscid flow vortex lines remain vortex lines, the 
above theorem is an statement about their evolution in time when the flow is 
isotropic and incompressible. The same is also true for vortex sheets which will 
follow from the discussion in section 2. 

Remark 1.4. The reason for presenting various upper and lower bounds is for 
their potential in applications. For instance, with regard to the above remark and 
under the same conditions, one can compute 7 2 (t) as the ratio of the enstrophy 
at time t to time t 0 . 

Proposition 1.2 will give us “tight” upper and lower bounds for | | p ds) 

in an obvious way. This can be used partially to get information about the mo- 
ments of C t . 
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Theorem 1.2. For an isotropic, incompressible flow, let a p {t) = (wf + w£ + 


ty|)/3. Then, for p > 0, 

(C t0 ) p < <(C t ) p > 

(0 

for 0 < p < 1, 


a p (t)[mm (£[0 ,i]{|C , (s;to)r 1 >] C to < <(C t ) p > < (C t ) p < c p 

(Cto) p (tt) 

for p > 1, 


(C,.)» < c? (C..)' < (C,)" < ((C.)') < 3o,(()(max, < |0, 1 |{|C'( 51 t 0 )|'’-‘}] C«„ 

(*0 


where C 2 (ci) is the minimum (maximum) of the coefficients of C to in the upper 
(lower) bounds for (C t ) in Theorem 1.1, and without loss of generality, we have 
assumed [o,6] = [0,1]. 


Proof. The easiest way to handle (») is to raise the left hand side of (t) in 
Theorem 1.1 to power p, take p inside the log , and then use (18) with <j> = exp(x). 
The left hand side of («) and the right hand side of (m) are also the cosequences 
of (18) with 4> = x p in an obvious way, and the remaining parts are easily followed 
by Theorem 1.1./// 

The next corollary will give us information about the moments of evolution 
of a straight line segment in turbulent flow and also the moments of dispersion 
of two points as time goes on. 

Corollary 1.2. Let di = mm[c5,Op(l)] , ^2 = max[c^,Za p {t)\ with ci,C 2 and 
a p (t) as in the above theorem.Let C(s;$o) = oi + s(a 2 — ai), se[0, l] . Then for 
an isotropic, incompressible turbulence, 

dx|o 2 - ai| p < ((C t ) p ) < d 2 |a 2 - Oi| p for p > 0, (t) 

(|z(o 2 ,f) - x(oi,f)| p ) < d 2 |a 2 - ai| p for p > 0. («) 

Proof. Note that | | = I 02 ~ <*i| and use the above theorem./// 

2. Surface Stretching 

In this section we will extend the above results to the evolution of a surface 
area in turbulent flow. It turns out, just as in the work of Cocke (1969), that 
only minor modifications are needed to do so. 

Let 6° l and 6°k be two infinitesimal material line at t — t 0 . We can form an 
infinitesimal material surface by taking the vector product of these vectors, i.e. 
6°S = 6°l X 6°k. This at time t becomes, 

3 3 

6S = Six 6k = (%2XjS 0 lj) x Xj6°kj). 

3 = 1 3 = 1 


( 14 ) 
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Let the matrices W, A and the eigenvalues, Wi’s, of W, be as before. From the 
identity, 

|t>i x v 2 \ 2 = |vi| 2 |« 2 | 2 - («i * w 2 ) 2 , (15) 

for any given vectors Vi and t >2 and the fact that |6°5| 2 = \A6°l x A6°k\ 2 , A 
being unitary, we can easily obtain, 

\6S\ 2 (A6°l x A6°k)l (A6°l x A6°k )% ( A6°l x A6°k)% 

|fi°S | 2 " W2Ws \AS°l x A6°fc| 2 + WlW 3 I A6°l x A6°k \ 2 + WlW 2 \A6°l x A6°k | 2 

= (sm 2 0cos 2 0)u>2t«3 + (sin 2 Osin 2 tp)wiWz + (cos 2 6)wyW 2 . 

( 16 ) 

Now (16) plays the role of (6), and we only need to modify Assumption (3) as 
follows; 

A ssumption; 

(4) wiw 2 , W 1 W 3 and w 2 w$ are identically distributed. ( Note that Assumptions 
(3) and ( 4 ) are equivalent for incompressible flows.) 

Now all we need to do is to replace wi,w 2 , by w 2 ws,w 5 wi,wiw 2 ; Sl,6°l 
by 6S,6°S , and C t ,C to by S t ,S to , and [a, 6] by D, respectively, in the above 
results, including the remarks, to obtain the new ones corresponding to surfaces. 
Where we let S(u, v; to) be a parametric representation of a nonrandom surface 
on a region D on the plain, S to be its area, and, 


'<-!L 


,aj(g(u,t>;t 0 ),t) dz(S(u,t;;to),t) , 


du 


dv 


dudv , 


(17) 


the area at time t. 

The only nonsymbolic modification of the proofs are: (a) In Theorem 1.2 the 
area of D has to be one or otherwise the right (left) hand side of (m) ((**)) has 
to be multiplied by that area to the power p — 1 due to the correct usage of 
Jensen’s inequality, (b) In Corollary 1.2 the notion of a distance between two 
points has to be replaced by an area of a region on a plane and its left hand side 
of (ii) to be interpreted correctly. 

Remark 2.1 We could have always used Proposition 1.2 and its counterpart 
for surfaces to obtain upper and lower bounds for moments of material lines, 
material surfaces, arc lengths, and surface areas at the expense of having different 
constants. Compare Proposition 1.1 with Proposition 1.2 when p = 1. The 
difference between these two propositions becomes significant if one can realize 
physically non-isotropic incompressible flows that satisfy only one pair of the 
assumptions, involved in these propositions, rather than all of them. Finally, 
for analogous results concerning homogeneous turbulence we invite the reader 
to consult Corrsin (1972). 
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Appendix 

(Jensen's Inequality) Let X be a random variable and (j> a convex (concave) 
function containing the range of X. Assume both X and 4>{X) have ensemble 
averages, then 

•KW) < (>) <*(*)>• (is) 

Proof. See any standard graduate textbook in probability theory or measure 
theory e.g. Billingsley (1986) p.283. 

The following special case of Jensen’s inequality has been used frequently; let 
Pi, P2, and P 3 be three positive numbers whose total sum is one, let aj., 02 , and 
03 be any real numbers. Then with <j> as above we have, 


3 3 

*£>«,) < (>) J^Pi^(a,). 

1=1 1=1 


(19) 


Proof. Let X in (18) be the random variable which takes the value a,- with 
probability p,-, t = 1, 2, 3./ / / 

(Dunford and Schwartz[ 5], p.535). Let (5,£,p) and (5i,Ej,pi) be positive 
measure spaces. Assume p(5) = 1 . Then if K is a p. x pi - measurable function 
defined on S x Si, 


/ exp{ I log\K(s,s 1 )\p(ds)}p 1 (dsi) < exp{ / \log{ / |/ir(s,si)|pi(dsi))]p(<fs)}. 

JSi Js Js J Si 

( 20 ) 
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Renormalization group analysis of turbulence 

By L. M. SMITH 


1. Objective 

The objective is to understand and extend a recent theory of turbulence based 
on dynamic renormalization group (RNG) techniques. The application of RNG 
methods to hydrodynamic turbulence has been explored most extensively by 
Yakhot and Orszag (1986). They calculate an eddy viscosity consistent with the 
Kolmogorov inertial range by systematic elimination of the small scales in the 
flow. Further, assumed smallness of the nonlinear terms in the redefined equa- 
tions for the large scales results in predictions for important flow constants such 
as the Kolmogorov constant. The authors emphasize that no adjustable param- 
eters are needed. The parameterization of the small scales in a self-consistent 
manner has important implications for sub-grid modeling. 

2. The RNG Transformation 

Renormalization group methods were first developed for quantum field theo- 
ries. They were later applied to the theory of critical points in materials that 
undergo phase transitions (Ma, 1976). Predictions for the universal exponents 
characterizing the behavior of thermodynamic quantities near critical points are 
quite accurate. The common feature of the physical phenomena amenable to 
RNG analysis is a lack of characteristic length and time scales. 

The lack of characteristic length and time scales in turbulence makes RNG 
methods attractive. The universality of the inertial range spectrum in widely 
varying turbulent flows is also suggestive. 

The RNG transformation consists of two steps. First, small scales are elimi- 
nated by an averaging procedure. Second, space is rescaled. New independent 
variables are defined on the original intervals by the rescaling. In most cases, 
the dependent variables must also be rescaled. 

A set of equations is renormalizable if it is unchanged by the RNG transfor- 
mation. Renormalizability implies scale invariance. Usually a set of equations is 
renormalizable only for specific values of its coefficients and the scaling parame- 
ters. These points are called fixed points. However, the physics of more general 
cases is often well described by the physics at a fixed point. 

The method of attack is to iterate the RNG transformation of the equations. 
With each transformation the coefficients in the equations change. One looks 
for a situation in which this iteration procedure converges. 

In addition to redefining coefficients of existing terms, the scale elimination 
often generates terms of different form than those in the original equations. 
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These new terms can be classified as irrelevant, marginal or relevant according 
to whether they decay, are constant or grow when rescaled. One can ask if a 
fixed point exists in the absence of new terms. If so, all new terms must be 
irrelevant for the system to be truly renormalizable at that point. 

3. The Basic Premise of the RNG Analysis of Turbulence 

The theory is based on the postulated equivalence between inertial range solu- 
tions of the Navier Stokes equations subject to initial and boundary conditions, 
and homogeneous isotropic flow driven by a Gaussian random force (Forster et. 
al., 1977, Yakhot and Orszag, 1986). The model equations are then 

^ + ( v • V)v = f - -VP + i/ 0 V 2 v (1) 

V • v = 0 (2) 

where v(x,t) is the velocity, P the pressure, p the density, v a the kinematic 
viscosity and f the forcing. The domain of equations (1) and (2) is unbounded. 

The white noise force is given by its correlation function in wavevector, fre- 
quency space. The correlation is assumed to obey a power law spectrum, 

(/<(k,w)/y(-k, -w)> a |k|~ v (3) 

where the brackets indicate an ensemble average. The exponent y is chosen to 
give the inertial range energy spectrum. Once y is fixed, there are no adjustable 
parameters in the problem. 

4. A Revised RNG Analysis 

Yakhot and Orszag show that analysis of (l)-(3) using the full RNG trans- 
formation yields the scaling laws of velocity correlations, and thus the energy 
spectrum. If y is set equal to the number of dimensions, 3, the Kolmogorov 
spectrum is recovered: E{k) oc k ~ 4 * 6 / 3 where k = |k|. Amplitudes, however, are 
left undetermined. 

By performing only the scale elimination, and abandoning the rescaling, they 
are able to find both scaling laws and amplitudes. Then E(k) = Ko,Rw G e 2 / 3 A: _6 / 3 , 
where e is the dissipation rate and Korng is the RNG prediction for the Kol- 
mogorov constant. 

Rescaling is used only to justify neglect of new terms generated by the elimi- 
nation procedure. The terms of concern are cubic in the velocity vector and are 
marginal with respect to the fixed point found in their absence. One wonders 
how the results would change if the cubic terms are retained. A goal of the 
present research is to assess the effect of these terms on the system. 
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5. The Effect of the Small Scales 

The theory developed by Yakhot and Orszag is an attempt to calculate the 
effect of the small scales on the large scales in turbulence. Their method deter- 
mines that the large scales ‘feel’ the small scales as an eddy viscosity. 

Equations (l) and (2) are written in wavevector, frequency space and the 
pressure is eliminated by taking the curl of the curl. The equations for the 
Fourier coefficients of the velocity field are then expanded in a power series via 
the introduction of an ordering parameter which multiplies the nonlinear term. 

A narrow band of wavenumbers is removed by averaging over their force field. 
The averaging procedure replaces the contribution of the nonlinear interaction 
of those wavenumbers with a term linear in the velocity vector for the remaining 
wavenumbers. The nonlinear interaction is only approximately represented in 
this term. The approximation is due to truncation of the power series at sec- 
ond order and neglect of terms cubic in the velocity vector for the remaining 
wavenumbers. 

The coefficient of the linear term is an integral. The integral is evaluated in 
the limit 0 ~ A; << k c , where k c is the low wavenumber cutoff of the eliminated 
band. In this limit the integral is proportional to k 2 . Thus the large scales see 
the small scales as (approximately) a viscous term. 

Iteration of the elimination procedure produces an equation for the large scales 
and long times identical in form to equation (l). The molecular viscosity u 0 is 
replaced by an eddy viscosity, 

_ \ 

— + A(v • V)v = f - -VP + z/ r V 2 v (4) 

ot p 

where l/t is the eddy viscosity and A is the nondimensionalized ordering param- 
eter. The eddy viscosity depends on k Ci 

"r = Vo[ 1 + P 3 t(k~* - K 4 )] 1/3 (5) 

where /? is a function of A, k c is the last eliminated wavenumber and k,i is the 
viscous cutoff. If one now takes the limit k c ~ k ~ 0, 

i>t ~ ( 6 ) 

In this limit the ‘renormalized’ equation (4) has a fixed point, at which all 
subsequent analysis takes place. 

6. Evaluation of Universal Flow Constants 

The nondimensionalized ordering parameter A is proportional to (y + l) 1 / 2 , 
where — y is the power of k in the force correlation function (see (3)). Recall 
V — 3 for a Kolmogorov inertial range. Using y = 3, the fixed point value of A 
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is 11.5. However, the RNG results depend on the simultaneous extrapolation of 
k c ~ k ~ 0 and y approaching -1. 

The assumption of small (y + l) allows a power series expansion of the renor- 
malized equations. The equation for the zeroth order velocity coefficient is the 
Langevin equation, 


(-tw + u T k 2 )\^ (k, w) = f , (7) 

where v(k,w) is the Fourier amplitude, v = v^°) + Av^ + A 2 v^ 2 ) -|- ..., k and 
u> are small, and l>t is given by equation (6). The value of 0 at the fixed point, 
evaluated to zeroth order in A, is found to be 0 = .493. 

Evaluation of the Kolmogorov constant follows a power series expansion at 
the fixed point of the energy equation (the equation for the two-point velocity 
correlations). The analysis is again for fc c ~ fc ~ 0 with Ur given by (6). 
When terms to lowest nontrivial order are retained, its value is found to be 
Korng — 1.617 (Dannevik et. al., 1987). 

The same scale elimination procedure can be performed on the equations for 
the advection of a passive scalar. The renormalized equations for the large 
scales are characterized by an eddy diffusivity. A power series expansion with 
fc c ~ fc ~ 0 leads to a prediction for the Batchelor constant, Ba. 

7. Goals 

a. Further steps in the analysis remain to be understood. For example, the 
RNG k — e model leads to a prediction for the von Karman constant. 

b. One would like to assess the importance of the neglected new terms gener- 
ated by the scale elimination procedure. These terms are cubic in the transform 
coefficients of the velocity vector for the large scales. 

c. An exploration of other forcing functions will be revealing. 

i. We have already investigated the possibility of a non-white noise, power 
law spectrum. The force correlation is allowed to fall off with frequency and 
considered proportional to k~ v u~ b . We find that there is no positive value of 6 
consistent with both dimensional analysis and a Kolmogorov energy spectrum. 

ii. A completely satisfactory theory of turbulence should account for spatial 
intermittency. In the present RNG formulation, intermittency might be included 
by changing the statistics assumed for the forcing function. 

d. It should be possible to extend the theory to flows other than homogeneous, 
isotropic turbulence. Perhaps one can do away with the artificial forcing for 
realizable kinds of turbulence. 

e. Finally, the RNG sub-grid and k — e models should be tested. 
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Low-dimensional chaos in turbulence 

By J. A. VASTANO 


Abstract 

Direct numerical simulations are being performed on two different fluid flows in 
an attempt to discover the mechanisms underlying the transition to turbulence 
in each. The first system is Taylor-Couette flow; the second, two-dimensional 
flow over an airfoil. Both flows exhibit a gradual transition to high-dimensional 
turbulence through low-dimensional chaos. The hope is that the instabilities 
leading to chaos will be easier to relate to physical processes in this case, and 
that the new understanding of these mechanisms can then be applied to a 
wider array of turbulent systems. 


In the past decade a new understanding of nonlinear processes in nature has 
been provided by the application of mathematical ideas on chaos and bifurcation 
theory (Ott 1981, Swinney 1983). An important result for turbulence researchers 
has been the finding that some systems undergo a transition to turbulence by 
first becoming low-dimensionally chaotic (Brandstater and Swinney 1987, Keefe 
1987, Sreenivasan 1985). Turbulence models should benefit from a better un- 
derstanding of weakly turbulent cases in which low-dimensional models can in 
principle capture all of the system dynamics. 

The purpose of this research project is to study instability mechanisms in 
two model cases that have exhibited a gradual transition to turbulence, i. e., 
a transition to turbulence via low-dimensional chaos. The hope is that the 
instabilities leading to chaos will be easier to relate to physical processes in 
these models, and that the new understanding of these mechanisms can then be 
applied to a wide array of turbulent systems. 

The first of the two projects, undertaken with the advice and assistance of 
Drs. Robert Moser (NAS A- Ames) and Lawrence Keefe (CTR), is the com- 
putation of the Lyapunov exponent spectrum for a model of Taylor-Couette 
flow, the flow between concentric rotating cylinders (Brandstater and Swinney 
1987, Moser 1983). Experiments on moderate aspect ratio Taylor-Couette sys- 
tems have observed a transition from laminar, time-independent Couette flow 
to periodic, quasiperiodic, and finally to low-dimensional chaotic motion. The 
physical mechanisms driving the low-order transitions have been found, but the 
transition to chaos is not well understood. The numerical method that we will 
use was developed by Moser (1983) to study curved channel flow and assumes 
axial periodicity. The method recovers the low-order transitions well, but the 
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axial constraint may significantly affect the chaotic states. The experimental 
evidence in the chaotic regime includes measurements of wavespeeds associated 
with dominant frequency components in the flow as well as estimates of the frac- 
tal dimension for the chaotic attractor. We can, therefore, verify the accuracy 
of our simulation: the wavespeeds cam be measured just as they are in experi- 
ment, and the fractal dimension can be computed directly from a knowledge of 
the Lyapunov exponent spectrum for the flow (Keefe 1987). If our numerical 
results differ greatly from the experimental results on wavespeed, we may need 
to consider a finite axial extent model. Given an good numerical model of the 
flow, our dimension estimate will be more accurate (as accurate as computer re- 
sources allow) than dimension estimates obtained indirectly from experimental 
time series data. 

Although the dimension estimate we obtain from the Lyapunov exponent spec- 
trum is necessary to validate the model (or, perhaps, the experiment), there is 
a great deal more to be learned from the exponents. The Lyapunov exponents 
for a system measure the exponential rates of growth or decay for infinitesimal 
perturbations in each phase space direction. There are, therefore, an infinite 
number of Lyapunov exponents for any spatially-extended system, but for dis- 
sipative systems almost all of the exponents are large and negative, signifying 
the decay of transients to some attractor. The number of positive Lyapunov 
exponents is a rough indication of the number of distinct mechanisms destabiliz- 
ing the flow. Moreover, the time-varying contributions to the long-time average 
exponents will indicate when in the system evolution the chaos-generating mech- 
anisms are operating. The spatial structure of the eigenvectors associated with 
the positive Lyapunov exponents at those times may indicate the nature of that 
mechanism by showing where in the flow it is operating. 

Progress on this work since 1 Sept 1988 includes familiarization with the NAS 
operating environment and the Vectoral programming language, and full imple- 
mentation of the curved channel code on the Cray 2. The additions to the code 
necessary to compute the Lyapunov exponents have been inserted and tested 
on periodic and quasiperiodic Taylor vortex flow. The basic code (without the 
exponent calculations) is being run in order to estimate the minimum resolution 
needed to model chaotic Taylor-Couette flow. Moser’s previous results showed, 
for example, that a 32 x 16 x 32 (R — d — Z) simulation of a periodic (wavy 
Taylor vortex) flow reproduced the azimuthal wave speed found in experiments 
to within the experimental uncertainty. Lowering the resolution to 16 x 16 x 16 
only changes the wave speed by 1%, while a resolution of 16 x 16 x 8 gives a 5% 
deviation. Moser attained good agreement with experiments on quasiperiodic 
(modulated wavy Taylor vortex) flows using a 64 x 32 x 64 model. We are cur- 
rently checking the effects of lowered resolution on this case. It seems likely that 
only moderately-high resolution will be needed for chaotic solutions at Reynolds 
numbers close to the quasiperiodic-chaotic transition. 

A second, related research project is being conducted in collaboration with 
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Dr. Thomas Pulliam (NAS A- Ames) (Pulliam 1989, Pulliam and Vastano). This 
is a study of the transition to turbulence in two-dimensional flow over airfoils 
at high angles of attack. Experiments on forced flow over airfoils have found 
that a transition to a chaotic state can be produced (Stuber and Gharib). Nu- 
merical simulations of two-dimensional unforced flow had found a transition for 
increasing Reynolds number from periodic flow to aperiodicity that included one 
period-doubling (Pulliam 1989). We have now seen a period-doubling cascade 
of bifurcations and a transition to low-dimensional chaos past a period-doubling 
accumulation point. We are currently attempting to characterize the observed 
chaos using Poincar^ sections and Lyapunov exponent estimates from time se- 
ries data. As in the Taylor-Couette research, the goal here is to identify the 
physical mechanisms causing the period-doubling to chaos, rather than simply 
establishing that such a transition exists. 
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Close interactions of 3D vortex tubes 

By M. V. MELANDER 1 

Abstract 

The motivation for studying close vortex interactions is briefly discussed in 
the light of turbulence and coherent structures. Particular attention is given to 
the interaction known as reconnection. Two reconnection mechanisms are dis- 
cussed. One is annihilation of vorticity by cross-diffusion, the other is an inviscid 
head-tail formation. At intermediate Reynolds numbers both mechanisms are 
operating. 

1. Introduction 

Close interaction and self-deformation of vortices are intriguing facets of fluid 
mechanics. These interactions occur repeatedly in turbulent flows. We can, 
therefore, expect that insight into fundamental vortex interactions will lead to 
a deeper understanding of turbulence. The theory of turbulence in 2D flows 
illustrates this point clearly. In 2D flows fundamental vortex interactions such 
as axisymmetrization, merger, and dipole formation provide the mechanisms 
for the emergence of large scale coherent structures (McWilliams 1984). These 
structures in turn invalidate earlier statistical theories based on random phases 
(Babiano et al 1987). Compared to 2D vortex dynamics, our understanding of 
three dimensional vortex interactions is in its infancy. In order to improve this 
level of understanding, it seems logical to begin by investigating the most basic 
vortex interactions. This report marks the beginning of a thorough investigation 
of such interactions. 

We confine the study to unforced incompressible flows with no density vari- 
ations and unbounded domains. Turbulence may then be viewed as a tangle 
of interacting vortices. This follows directly from the integral formulation of 
Navier-Stokes equations. 


1 f w(r') x (r - r ; ) , 

4 J |r-rf V + V * 

(1) 

u) = V x u 

(2) 

Du) 

— = (w • V)u + I/Aw 

(3) 


1 Permanent address: SMU, Texas 
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FIGURE 1. Schematic for the calculation of the velocity field u(r) generated 
by a vortical structure. 

When the fluid is at rest at infinity, the potential component in (l) vanishes. 
The Biot-Savart law (l) then shows that the vorticity governs the flow. The 
vorticity transport equation (3) in turn describes how the vorticity evolves with 
the flow. Some attractive features of the formulation (l) - (3) are: I) The 
irrotational part of the fluid does not appear in this description. This fact is the 
motivation behind computational methods known as “vortex methods” (Leonard 
1985). Unfortunately, these methods make strong assumptions concerning the 
local vorticity distribution and may need further development before they can 
safely be used to investigate delicate vortex interactions. II) The description (l) - 
(3) appeals to our physical intuition. It is conceptually easy to give a qualitative 
prediction of the short time evolution of a given vorticity configuration by means 
of arm-waving and finger-twisting arguments. It is also possible to understand 
vortex dynamics in terms of primitive variables and force-balancing (Moore and 
Saffman 1972); although this approach leads to the same physical results, it 
seems considerably more complicated. Ill) A systematic description of how the 
small scale vortical structures influence the far-field u(r) can be obtained from 
the Biot-Savart law (l). 

Let us consider the last point in detail. Let a vorticity configuration be given 
near position r' (Figure 1). The Biot-Savart integral gives the velocity field 
u(r) generated by this vorticity configuration at position r. The integration 
extends over the entire vorticity configuration. Clearly the integrand depends 
on r, so it is an enormous task to evaluate the integral for many different r 
positions. However, if jr' — r| is large, we can expand the integrand in inverse 
powers of |r' — r| . The coefficients in this expansion are moments of the vorticity 
configuration and hence independent of r. These moments describe the shape 
and internal structure of the vorticity configuration. Far away, only the lowest 
order moment matters, but as we let r approach the vorticity configuration, an 
increasing number of moments must be included. We see that physical space 
moments yield a systematic description of internal structures influence on the 
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far-field. This idea has been used with success in 2D (Melander et al 1986). 

The above discussion shows that the interaction between well separated vor- 
tices can be calculated efficiently in a straightforward manner. However, when 
jr — r' | is small, we are faced with an entirely different situation. We no longer 
have a natural expansion parameter for the evaluation of the Biot-Savart in- 
tegral. The integrable singularity in (l) tends to amplify the influence of the 
local structure. As a consequence, the local self-induced velocity can become 
very large for slender vortices. Such slender vortices can easily develop through 
vorticity stretching. A more subtle effect of small scale vorticity lies in its ability 
to alter neighboring large scale vortical structures. A concrete example from 2D 
is the axisymmetrization of a single isolated vortex (Melander et al 1987); here 
the small scale structures (filaments) influence the large scales in such a way as 
to make the vortex core circular symmetric. Because of these effects, we must 
treat the local vorticity structure and the associated small scales carefully. A 
thorough study of close vortex interactions via direct numerical simulations will, 
therefore, be most helpful for modelling of the local vorticity distributions and 
the associated small scales. 

Vortex dynamics reveals its most useful role when viewed in the light of co- 
herent structures. These vortical structures have been observed in many flows 
that were previously regarded as fully random (Hussain 1986) . Often these large 
scale structures are superimposed by fluctuating vorticity (incoherent vorticity), 
which makes the structures less obvious, but not less important. In spite of the 
incoherent vorticity, the large structures still dominate the far-field. The inco- 
herent vorticity does not influence the far-field directly, but can have an indirect 
influence by altering neighboring structures, as discussed above. Many inter- 
esting questions arise concerning the formation, persistence, characterization, 
topology, and interaction of coherent structures. Some of these questions are 
best approached through vortex dynamics. In this report we concentrate on in- 
teractions that can change the topology. These interactions are generally known 
as reconnections, although other names such as cut-and-connect, cross-linking, 
and fusion are also used. We examine a number of reconnection interactions 
under very idealized conditions. The interactions are simulated in isolation such 
that non-local effects do not obscure the dynamics. Also, there is no incoherent 
vorticity in our initial conditions. 

2. Initial Conditions, Simulation Methods, and Diagnostics 

The simulations were performed using a dealiased spectral (Galerkin) method 
with a fourth order predictor-corrector algorithm for the time advancement. 
This algorithm solves the Navier-Stokes equations (or the corresponding equa- 
tions with superviscosity) in a cube with periodic boundary conditions. The 
evolution of a passive scalar at a unity Schmidt number is calculated simultane- 
ously. The initial conditions (Figure 2a-f) were chosen as rectilinear vortices with 
a Gaussian vorticity profile. The evolution starting from the initial condition 
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strong 



FIGURE 2. Initial conditions for the simulations. 
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shown in Figure 2a was simulated using both Newtonian viscosity (Re = 2000) 
and superviscosity, while in the cases (b-f) only superviscosity was used. The 
spatial resolution in all these cases was 96 meshpoints in each direction. The 
simulation starting from the initial condition shown in Figure 2g is discussed at 
great length elsewhere (Melander and Hussain 1988 a,b,c). 

The simulations produced databases which were analyzed on graphics worksta- 
tions using the “interactive” program TURB3D. Displays of the spatial vorticity 
distribution at various stages of the evolution were crucial in obtaining a physical 
understanding of the interactions. Cross sections of vorticity and velocity fields 
were also very helpful, whereas displays of other quantities such as helicity, dis- 
sipation, and enstrophy production were found to be of secondary importance. 
Tracing of vortex-lines turned out to be a major disappointment, for it seemed 
impossible to find just the right set of vortex lines to clearly illustrate a given 
interaction. In general, a few vortex-lines were not enough to reveal the dynam- 
ics, while a larger number such as 20 gave a picture that was too complicated to 
follow. Perhaps displays of vortex tubes via Clebsch potentials would be more 
appropriate. 

3. Discussion 

In order to achieve a higher effective Reynolds number, a superviscosity was 
used in most simulations. The superviscosity leaves the largest scales more in- 
viscid than does a corresponding Newtonian viscosity. However, it is also known 
that the superviscosity generates artificial small scale structures. Comparison of 
two simulations, starting from the same initial condition (Figure 2a) but with 
different diffusion operators, showed similar evolutions. The most important 
difference between these two simulations was the evolution time-scale; it was 
shorter in the superviscosity simulation. A similar effect has also been observed 
by decreasing the Newtonian viscosity, and hence increasing Reynolds number. 
We therefore conclude that the use of a superviscosity mimics a higher Reynolds 
number flow and that the artificial small scale structures do not significantly 
influence the large scale structures. 

Inviscid filament calculations with regularized circular cores show that two 
orthogonally offset vortices become locally antiparallel (Schwarz 1983 and 1985, 
Siggia 1985). Careful studies show that if the same filament calculations are 
continued further in time, then a finite-time singularity develops (Pumir and 
Siggia 1987). This singularity is most likely caused by insufficient degrees of 
freedom in the vortex cores, that is, inadequate modelling. Nevertheless, some 
researchers feel that a similar singularity might be present in the full Euler equa- 
tions (Pumir and Kerr 1987, Kerr and Hussain 1988). It has also been suggested 
that this singularity is present in the Navier-Stokes equations at sufficiently high 
Reynolds number (Pumir and Siggia 1987). Direct numerical simulations have 
failed to capture signs of this singularity. 
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FIGURE 3. Evolution of the initial condition shown in Figure 2a at Reynolds 
number 2000. The panels show isovorticity surfaces at 56% of the initial peak 
value. Time is measured in units of w max (0)/20. 
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Since it is almost impossible to prove or disprove the presence of a singularity 
numerically, a more fruitful avenue is to concentrate on understanding the op- 
erating physical mechanisms. One can argue that if it is possible to explain the 
evolution without including new terms in the governing equations, then there 
probably is no singularity. This argument rests on the general experience that 
a singularity is a sign of missing physics in the governing equations. 

The evolution starting from the initial condition shown in Figure 2a highlights 
the soft and pliable nature of vortex cores, thereby clearly showing that vortex 
cores do not remain nearly circular. The close proximity of the vortices results 
in pulling of hairpins from the outer layer of the vortex tubes (Figure 3). These 
hairpins form even before the main vortex tubes become locally antiparallel. 
Later more hairpins form, especially in the wake of the locally antiparallel vortex 
pair (hereafter referred to as the dipole). 

The dipole is initially squashed by self-induction (due to the curvature of the 
antiparallel vortex pair) and thereby form a head-tail structure. The dipole 
is propagating by mutual induction. However, an oppositely directed velocity 
field is generated by the remaining part of the vortex tubes. At first the dipole 
velocity is the largest, but later the dipole-strength diminishes and as a result 
the dipole is washed backwards. Thereby the curvature of the locally antiparallel 
vortices reverses, the self-induced motion, therefore, also reverses, and the dipole 
separates slowly. This reversal is the crucial part of the reconnection process 
and occurs at the same time as the new “reconnected” vortices begin to separate 
(Figure 3). 

The reversal is caused by the diminishing dipole-strength, which can happen 
by both a diffusive and an inviscid mechanism. The diffusive mechanism is 
annihilation of antiparallel vorticity by cross-diffusion in the contact-zone and 
is clearly important for low to medium Reynolds numbers. The result is a 
true reconnection, albeit this reconnection is only partial as the cross-diffusion 
is arrested by the reversal of the dipole propagation direction (Melander and 
Hussain 1988a). The inviscid mechanism is the head-tail formation. The dipole 
propagation velocity is determined almost exclusively by the head; therefore, the 
head-tail formation effectively diminishes the circulation in the dipole. In the 
high Reynolds number limit, this may result in an apparent reconnection. That 
is, the large scale vortical structures appear as though a topological reconnection 
has occurred, where in fact only a topology preserving entanglement has taken 
place. The reasons for this conjecture are as follows. 

Several numerical simulations with different Reynolds numbers and different 
diffusion operators show that the large scale vortical structures rearrange on 
a time-scale T, which does not increase with Reynolds number. Weakening of 
the dipole is responsible for this rearrangement. If viscous effects are solely re- 
sponsible for this weakening, then \uAu\ must be at least finite as the Reynolds 
number tends to infinity. This in turn implies that | Aw| must become unbounded 
within the finite time T, which can only happen through axial stretching of the 
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FIGURE 4. Asymmetric entanglement and apparent reconnections; the cir- 
culation ratio is four. The isosurface level and timescaling are the same as in 
Figure 3. 
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FIGURE 5. Sketch of Husain and Hussain’s rectangular jet experiment. 

dipole. This stretching is caused by lengthening of the antiparallel vortices in 
the contact-zone and by the strain form outside the contact-zone. From axisym- 
metric vortex dynamics (Stanaway et al 1988), we know that the lengthening 
effect is insufficient to cause a finite time blow up of |Aw|. The external strain 
is also bounded due to the finite length of the contact-zone. 

Neither of the simulations starting from the initial condition shown in Figure 
2a had sufficient resolution to clearly separate the inviscid and viscous mecha- 
nisms. It does not seem likely that sufficient resolution can be obtained in the 
near future. In this respect the asymmetric initial condition (Figure 2b) is more 
promising due to the fact that the circulation ratio can be made small such as 
to allow for combined simulation and analysis. The dynamics of this simulation 
is also much easier to grasp, see Figure 4. The weak vortex wraps around the 
stronger one and thereby forms a large hairpin. Vortex stretching occurs mainly 
in the legs of the hairpin. The tip of the hairpin is unstretched and dispersed 
in the azimuthal direction. Therefore, the curvature of the legs accounts for 
most of the self-induced motion, which is such as to split the hairpin into two 
diverging helical structures. 

The evolutions of the initial conditions shown in Figures 2c-f were also in- 
vestigated and explained. I later became aware that the latter occurs in a 
rectangular jet experiment by F. Hussain, see Figure 5. Pulsing of the jet gen- 
erates rectangular vortex rings, which initially undergo a near-recurrence with 
90 ° axis-flip. This phenomenon can be explained by a filament model (Bridges 
and Hussain 1988), except for rapid mixing processes near the corners. Further 
downstream an overtaking collision occurs and the resulting reconnections pro- 
duce four rings. The overtaking collision occurs because the rectangular vortex 
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FIGURE 6. Simulation of an overtaking collision; the fast vortex pair has 
twice as much circulation as the slow pair. The isosurface level and timescaling 
are the same as in Figure 3. 
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rings do not comply with the traditional leap-frog process of axisymmetric vor- 
tex rings. The experiment compares well with the simulation shown in Figure 6 
and clearly validates our selection of idealized initial conditions. 

4. Conclusion 

The insight gained by viewing the evolutions of the above idealized initial 
conditions on graphics workstations has proved to be invaluable in the process 
of constructing mathematical models of the evolutions. The modelling, which is 
now under way, centers around two types of initial conditions. One type is the 
symmetric configuration of two antiparallel vortices with sinusoidal perturba- 
tions (Figure 2g). The analysis of this model also applies to the collisions shown 
in Figure 2e-f. The asymmetric entanglement shown in Figure 2b and Figure 4 
is also being modelled for small values of the circulation ratio. 
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Numerical simulation of viscous vortex rings 

By S. K. STANAWAY 

1. Motivation & Objectives 

This work is directed toward understanding vortex interactions and their role 
in turbulent flow. The objectives are twofold. First, to use the existing ax- 
isymmetric code to study the annihilation process of colliding vortex rings and 
determine the relevance of this problem to similar 3D phenomena. The second 
objective is to extend the code to three dimensions. The code under development 
is unique in that it can compute flows in a truly infinite domain (i.e. without 
periodic boundary conditions or approximations from truncating the domain). 
Because of this, we are able to compute the feu: field sound, and therefore, con- 
tribute to improved models of turbulence generated noise for this class of flows. 
Issues which can be addressed by the code include: effects of viscosity on mode 
selection in azimuthal breakdown of vortex rings (i.e. the Widnall instability); 
reconnection, the associated production of small scales, and the time scale of the 
process. 

2. Accomplishments 

Objective 1: Study head-on collision of viscous vortex rings 

During the CTR Summer Program, a collaboration was undertaken between 
the author, K. Shariff and F. Hussain to study the head-on collision of two 
identical axisymmetric viscous vortex rings through direct simulations of the 
incompressible Navier-Stokes equations. The result of this effort is described 
in the proceedings for the 1988 Summer Program of the CTR (Stanaway et al. 
1988). 

Initial conditions of varying core shapes (thin rings, and Hill’s spherical vortex 
rings) and Reynolds numbers (350 to 1000) were considered and the subsequent 
annihilation process was studied through time history of circulation, vorticity 
contours and local contributions to annihilation. The results provide a database 
to model such a problem and also information about the level of detail required 
to accurately model this phenomenon. A typical interaction proceeds as follows. 
Two vortices of opposite sense and the same strength approach each other by self- 
induction, the radii increase by mutual induction, and vorticity cancels through 
viscous diffusion across the collision plane. Following contact, we observed (for 
the cases considered here) that the vorticity distribution in the core forms a 
head-tail structure, a behavior which has also been seen in inviscid calculations 
(Shariff et al. 1988), 3D viscous calculations with periodic boundary conditions 
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(Melander & Hussain 1988, MH), and experiments (Oshima 1978). In examining 
the local contributions to the total vorticity annihilation in both the head and tail 
regions, we conclude that the tail does contain a significant part of the vorticity 
and should be included in an accurate model of the annihilation process. 

The characteristic time of vorticity annihilation is compared with that of a 3D 
collision experiment (Schatzle 1987) and 3D numerical simulations (MH). It is 
found that the annihilation timescale for the axisymmetric collision is faster than 
the viscous timescale, a 2 /i/, and slower than the timescale set by the circulation, 
°o/(ro*') 1/2 , where a Q is the initial core radius, r o is the initial circulation, 
and v is the viscosity. This indicates that the local effects are important in 
enhancing annihilation, however, nonlocal effects such as vorticity realignment 
are also important in 3D. One might expect that during the initial stages of 
the collision, local effects are dominant, and as the circulation in the symmetry 
plane weakens, the bridges strengthen and the out-of-plane strain becomes the 
more important effect. 

The flow is also computed to the large time Stokes flow limit where the circu- 
lation decays as t -3 / 2 and the vorticity distribution agrees with the quadrupole 
solution of the Stokes equations. In this limit, the self-annihilation is exactly 
twice the mutual annihilation. For one of the cases computed, the far-fxeld 
quadrupole sound is compared with the experimental results of Kambe & Minota 
(1983). The agreement is quite good even though the Reynolds numbers are very 
different. 

Objective 2: Extend axisymmetric spectral code to three dimensions 

An operational axisymmetric code which solves the Navier-Stokes equations 
in an unbounded domain is being modified to compute 3D flows. The method 
(Stanaway, Cantwell & Spalart 1988) uses divergence free basis functions, hav- 
ing the advantages of reducing the number of unknowns from four to two, and 
eliminating the need to solve a Poisson equation for the pressure at each time 
step. The solution is expanded in polar coordinates with the associated Leg- 
endre functions in the polar direction and Jacobi polynomials matched to an 
algebraic mapping of the radial coordinate. The third dimension, the azimuthal 
direction, is expanded in Fourier series. The above functions are used in order 
that the resulting matrix equations are numerically attractive; they are com- 
pletely orthogonal in two directions and have a constant and relatively narrow 
bandwidth in the third direction (bandwidth < 11). In addition, in the present 
method the matrices are symmetric and positive definite. 

The unknown coefficients in wave space, referred to as the + and — modes 
(in keeping with the notation of similar approaches), each have an ordinary dif- 
ferential equation describing their evolution. The axisymmetric basis functions 
are a subset of the three dimensional functions, meaning that in extending the 
axisymmetric code to 3D, much of the code remains intact. Specifically, the 
evolution equation for the + modes is unchanged. In the 3D case, an additional 
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evolution equation is required for the — modes. In the development of the ax- 
isymmetric code, it was necessary to compute the matrix elements analytically 
rather than numerically in order to minimize errors leading to ill conditioned 
matrices. This process, which was rather involved due to the complex coordi- 
nate system and mapping, was made possible with the aid of MACSYMA, an 
algebraic manipulation program. It is worth noting that the effort required to 
compute the elements of the matrices for the — modes was considerably less than 
it was for the + modes. The major effort, therefore, in extending the code to 
three dimensions has been in developing the transforms to and from wave space 
for the 3D basis functions. This step is nearly finished and will be tested in a 
straightforward manner by starting with an initial vorticity field, transforming 
to wave space, and then back to real space. 

So the problem of extending the code involves four phases: 

(i) Deriving the basis functions and evolution equations. 

(ii) Forming the matrices for the — mode evolution equation. 

(iii) Developing and testing three dimensional transforms going to and from 
wave space. 

(iv) Running test cases. 

Considerable progress has been made in extending the code to three dimen- 
sions and it is expected that the 3D code will soon be operational. 

3. Future Plans 

This new method gives us the unique capability to study many important and 
often controversial phenomena accurately, notably 

• the Widnall instability of an almost axisymmetric ring, in particular, the mode 
selection process; 

• viscous connection and splitting; 

• redistribution of vorticity by colliding or pairing vortex rings; and 

• generation, intensification and annihilation of vorticity through vortex ring 
interactions. 
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The effects of particle loading on 
turbulence structure and modelling 

By K. D. SQUIRES AND J. K. EATON 

1. Introduction 

Conventional wisdom implies that the presence of particles provides an ad- 
ditional dissipation or attenuation of turbulence. However, it is not clear how 
this extra dissipation may be incorporated into turbulence models or how it de- 
pends on the parameters of the problem including the particle time constant, 
the mass loading, and the various dimensionless parameters describing the tur- 
bulence. Simply focusing on the extra dissipation completely neglects the effect 
of turbulence structure on the instantaneous particle concentration field and 
the possibility of interactions between the particle cloud and instability mecha- 
nisms generating turbulence. For example, particles selectively concentrated in 
particular structures may cause rapid attenuation of that structure or trigger a 
new instability mechanism. Therefore, it is important to determine the effect of 
turbulence structure on the behavior of the particle cloud. 

The objective of the present research was to extend the DNS approach to 
particle-laden turbulent flows using a simple model of particle/flow interaction. 
The program addressed the simplest type of flow, homogeneous, isotropic tur- 
bulence, and examined interactions between the particles and gas phase tur- 
bulence. The specific range of problems examined include those in which the 
particle is much smaller than the smallest length scales of the turbulence yet 
heavy enough to slip relative to the flow. The particle mass loading is large 
enough to have a significant impact on the turbulence, while the volume load- 
ing was small enough such that particle-particle interactions could be neglected. 
Therefore, these simulations are relevant to practical problems involving small, 
dense particles conveyed by turbulent gas flows at moderate loadings. 

This report presents a sample of the results illustrating modifications of the 
particle concentration field caused by the turbulence structure and also illus- 
trating attenuation of turbulence by the particle cloud. 

2. Overview of the simulations 

The numerical method used in the present research is based on the pseudo- 
spectral technique developed by Rogallo (1981) for solving the incompressible 
Navier-Stokes equations. A modification of Rogallo’s code was used to simulta- 
neously track a large number of particles using a vectorized, second-order inter- 
polation scheme. This code is time advanced using second-order Runge-Kutta, 
and at each substep the fluid velocity at the particle position was calculated 
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using trilinear interpolation. Numerical experiments showed that more accurate 
interpolation schemes do not significantly improve the results. The new position 
was obtained by then integrating the particle equation of motion 


dvj(t) _ Uj(Xj(t),t) - Vi (t) 
dt t p 


( 2 . 1 ) 


dXjjt) 

dt 


= Vi(t) 


( 2 . 2 ) 


where X 3 (t) and v,(t) are the position and velocity of the particle, r p is the 
particle time constant, and Ui(X 3 (t),t) is the velocity of the fluid. Stokes’ law has 
been used to calculate the drag on a particle, and it is, therefore, assumed that 
the particle is much smaller than any lengthscale of the flow (the Kolmogorov 
scale, rj). It is also assumed that the particles are sufficiently dense such that 
other forces in the full equation of motion, e.g. buoyancy and added mass, are 
negligible compared to the Stokes drag. Finally, the particles are assumed to 
occupy a negligible volume fraction, and particle-particle interactions are also 
assumed negligible. 

To calculate the effect of the particles on the turbulence, the momentum 
equation of the fluid was modified by a source term on the right-hand side 


dui dui I dP d 7 Ui c . 

~xr + u i~aT ~ + — Z («» “ 

dt at poxi dxjdxj ■ r p p 

where c(x,-, t) is the particle mass per unit of volume, or particle density (repeated 
indices imply summation). 

The problem studied is homogeneous, isotropic turbulence in a box. “Natural” 
isotropic turbulence decays, and the lack of stationarity complicates analysis of 
the results. Therefore, the low wave number modes (large scales) are artificially 
forced to maintain stationarity using the scheme developed by Hunt, et al (1987). 
Results were obtained using 373000 particles for simulations using 32 s points and 
l.xlO 6 particles for simulations using 64* points. The relatively coarse grid 32 s 
calculations permitted a parameter study in particle time constant and mass 
loading ratio. For the 32* simulations the values of the particle time constant 
divided by the eddy turnover time used were 0.14, 0.75, and 1.50. The mass- 
loadings used in the 32* simulations were 0.1, 0.5, and 1.0. All three values of 
the mass loading were used for each of the time constants for these simulations. 
For the 64* simulations the values of the dimensionless time constant were 0.075, 
0.15, and 0.52 and the mass loadings were 0.1, 0.5, and 1.0. 

For each simulation the two phases were uncoupled and time-advanced to reach 
a statistically stationary state and eliminate initial transients of the forcing. At 
this point in time, the momentum equation of the fluid was coupled to the 
particle momentum equation with a specified mass loading. Once the phases 
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were coupled, the simulations were time advanced just over 12 eddy turnover 
times for the 32 s simulations and slightly over 6 eddy turnover times for the 64 s 
simulations. 

Results from the simulations of no coupling (4> = 0 where <j> is the mass loading) 
between the phases provide a baseline case to compare the effect of increased 
mass loading on the turbulence statistics. These simulations are also useful 
for examining the effect of turbulence “structure” on the particle concentration 
field. For these baseline simulations the Reynolds number based on the Taylor 
microscale was approximately 37. This relatively low value for the Reynolds 
number assured good resolution of the velocity field. 

3. Results and discussion 

Maxey (1987) has shown that particles will tend to accumulate in regions of 
high strain rate or low vorticity. This effect can be measured by using the second 
invariant of the deformation tensor, dui/dxj. The second invariant of dui/dxj 
is given by 

I Id = — ~(5 a — — Wj-Wy) (3-1) 

where s 1 is the magnitude of the strain rate, , u»y is the jth component of 
the vorticity vector and wyu/y is by definition the enstrophy. 

Thus, from equation (3.1) highly vortical regions correspond to I Id being large 
and positive (where the number density is small), while regions of high strain 
correspond to I Id being large and negative (where the number density is large). 

Contours of number density for each of the particles used in the simulations 
showed distinct regions where particles had collected. It was found that the 
peak number density is as much as 40-50 times the mean value for the lightest 
particles used in the simulations. By comparing these contours with the contours 
of the second invariant of dui/dxj , it was found that the particles showed a 
tendency to accumulate in regions of low vorticity or high strain rate. As the 
particles were made more sluggish (i.e., larger r p ), there is less of a tendency 
for them to accumulate in these regions. For the heaviest particle used in the 
32* simulations (r p /T e = 1.50), the particle concentration field was found to be 
essentially random. 

The preceding results were quantified by examining the conditional expecta- 
tion of the number density given the value of the enstrophy. The conditional 
expectation showed that the effect of increasing particle inertia reduced the like- 
lihood of finding the particles in regions of low vorticity (figure 1). For the 
64 s simulations, the lightest particles (r p /T e = 0.075) actually showed less of 
a tendency to accumulate in regions of low vorticity than did particles which 
were twice as heavy ( r p /T e — 0.150). This can be explained by the fact that 
very light particles can follow nearly all of the turbulent motions and would, 
therefore, show no preference to be in regions of low vorticity. 
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<p — 0.0 



FIGURE 1. Conditional expectation of number density given enstrophy for 
the 64* simulations, <j> = 0. 


For increasing values of the mass loading, the conditional expectation of num- 
ber density given the enstrophy for the light particles showed that they had less 
of a tendency to be in regions of low vorticity. However, this was not the case 
for the heavier particles. 

Correlations between the number density and enstrophy showed that as the 
loading was increased the lightest particles became less correlated with enstro- 
phy. For the larger time constants, however, the number density and enstrophy 
become more correlated for increasing values of the time constant and increased 
loading. Correlations between number density and II<x showed that as the mass 
loading was increased the correlation between number density and lid increased. 
This increase in the correlation was more significant for larger values of the time 
constant. 

Once the two phases were coupled, the turbulence evolved to a new stationary 
state within about two eddy turnover times. For all the time constants used in 
the simulation, it was found that the time required to come to a new equilibrium 
was longer for increasing loadings, and the maximum development time for any 
case was just over two eddy turnover times. Time averaged statistics such as 
turbulence energy showed that increased mass loading decreased the turbulence 
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FIGURE 2. Turbulence energy versus loading from the 32 s simulations, 
r p /T e = 0.75. 

energy (figure 2). These results are consistent with the fact that the drag of the 
particles on the turbulence acts as an additional source of dissipation. This can 
be shown be deriving the transport equations for (beginning with equation 
(2.3)). 

Frequency spectra measured along the particle path were nearly identical for 
all simulations. The differences in the turbulence spectrum for the three cases 
can be attributed to the fact that the lighter particles are found more often in 
regions of low vorticity and high strain rate than are the heavier particles. This 
causes a sampling bias resulting in the differences between the spectra along the 
particle path. Measured values of the particle velocity spectra were found to be 
in excellent agreement with a theoretical prediction of Csanady (1963). It was 
somewhat surprising to find such close agreement with Csanady’s theory in view 
of the fact that the particles are concentrated within particular regions of the 
turbulence. Examination of the power spectra of the turbulence for increasing 
values of the mass loading showed that energy was removed nearly equally from 
all frequencies. This may be due in part to the low Reynolds number of the 
simulations. Selective energy loss within certain frequency bands may occur at 
higher Reynolds numbers. 
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Correlations between enstrophy and pressure were decreased more by the light 
particles than by the heavy particles for increasing values of the mass loading. 
Since regions of low pressure are associated with regions of high vorticity, these 
results indicate that the lighter particles cause more distortion of these eddies 
than do the heavier particles. The fact that the lighter particles show a more 
pronounced tendency to accumulate in regions of lower vorticity (for tfi = 0) than 
do the heavier particles may provide some explanation for this. 

Correlations between I Id and pressure showed that for the time constants and 
mass loadings used in the simulations the correlation between these two quanti- 
ties remained reasonably constant, decreasing slightly for the lightest particles. 
Correlating lid with pressure correlates regions of low pressure with vortices 
(large positive lid) and regions of high pressure with straining regions (large 
negative I Id)- Therefore, these results indicate that the correlation between the 
straining regions with regions of high pressure must be increasing to compen- 
sate for the loss of correlation between the vortical regions and regions of low 
pressure for the lighter particles. 

In summary, it was found that for the case of zero loading there are significant 
effects of the turbulence “structure” on the resulting concentration field. These 
results were quantified by measuring conditional expectations and correlations. 
It is shown that the lighter particles show a strong tendency to be in regions of 
low vorticity and high strain rate. 

For increasing values of the mass loading the power spectra of the turbulent 
fluctuations do not show any frequencies to be preferentially modified by the 
particles. This may be partially due to the fact that there is not a large range 
of scales in the simulations. 

The lightest particles used in these simulations were found to modify the 
turbulence field differently than did the heavier particles. Since the light particles 
show a more pronounced tendency to accumulate in regions of low vorticity for 
<j> = 0, these regions are modified more by the light particles as the loading is 
increased than the heavy particles. Evidently, this selective modification by the 
light particles causes more of a distortion of the turbulent eddies than the more 
uniform modification by the heavier particles. 
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Analyses of turbulence structures in shear flows 

By M. J. LEE 

Outline of the research program and a recent progress in the studies of sheared 
turbulence are described. The research program reported here is directed at 
two goals: (i) understanding of fundamental physics of organized structures 
in turbulent shear flows; and (ii) development of phenomenological models of 
turbulence based on physical arguments. Three projects that have been carried 
out are: 

A. structure of sheared turbulence near a plane boundary; 

B. distortion of turbulence by axisymmetric strain and dilatation; 

C. study of energy transfer in turbulent shear flow. 

1. Motivation and Objectives 

Project A was motivated by the demonstration of remarkable similarity be- 
tween homogeneous turbulent shear flow and turbulent channel flow in instan- 
taneous velocity (and vorticity) fields and statistical correlations at the same 
dimensionless shear rate 5* = Sq 2 /e (Lee, Kim & Moin 1987). Here, S = dU/dy 
is the mean ‘shear rate,’ q 2 = uiu7 is twice the turbulent kinetic energy and e is 
the dissipation rate of the kinetic energy. The similarity in instantaneous vector 
flow topology as well as in statistical measures bears a profound implication 
towards a possibility of developing a ‘universal’ turbulence model. In particu- 
lar, this study showed that shear rate alone (without a wall) can produce the 
streaky structures, providing a strong support for the importance of shear rate 
in determining turbulence structures. 

The presence of a solid boundary affects turbulence in two fundamental ways: 

(i) generation of mean velocity gradient (via the no-slip condition) which, upon 
interaction with turbulence, supplies energy to it; and 

(ii) suppression of velocity fluctuations in its vicinity. 

It is hypothesized that these two effects are distinct in affecting turbulence struc- 
ture so that they can be accounted for separately. Project A centers on the anal- 
ysis of these two functions of a solid boundary and consists of three subprojects: 
uniform-shear boundary layer (USBL), plane Couette flow (PCF) and shear-free 
channel flow (SFCF). See figure 1 for schematic of the flows. 

The uniform-shear boundary layer (USBL) project aims at a theoretical ex- 
planation of the generation mechanism of the streaks in turbulent shear flows 
(Lee <k Hunt 1988). In USBL, the mean flow has a velocity profile U = Sy + Uo 
with the shear rate S uniform and constant, and the flow field is assumed to be 
inviscid with a slip velocity condition at the surface y = 0. The analysis was 
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FIGURE l. Schematic diagram of turbulent shear flows: (a) a uniform-shear 
boundary layer (USBL); (6) plane Couette flow (PCF); (c) shear-free channel 
flow (SFCF). 
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carried out for the evolution and distribution of turbulence by using rapid dis- 
tortion theory (RDT). Particular emphasis is placed on the differences between 
the shear-free boundary layer (see Uzkan & Reynolds 1967; Hunt Sc Graham 
1978; Hunt 1984) in redistributing turbulence near the surface by the blocking 
effect. 

The plane Couette flow (PCF) and shear-free channel flow (SFCF) projects 
were motivated by a need to generate numerical databases of turbulent wall- 
bounded flows with a constant total shear stress r(t/) = pdU/dy — puv = r w 
across the boundary layer (PCF) and zero mean velocity gradient near a wall 
(SFCF), respectively. These two cases are considered as building-block flows 
since most turbulence models are based on the wall-layer similarity with a con- 
stant shear stress r (Townsend 1976, §§5.4, 5. 7-5.9). It is of our interest to 
characterize organized turbulence structures in the near-wall regions at different 
mean shear rates and to examine the capability of the current scaling laws (or 
turbulence models) in describing turbulence characteristics of such flows (cf. El 
Telbany Sc Reynolds 1980, 1981). 

The purpose of Project B is to investigate the behavior of homogeneous tur- 
bulence subjected to axisymmetric distortion (contraction and expansion) and 
dilatation (Lee 1988). Closed-form solutions obtained by using RDT are com- 
pared with the numerical simulations (Lee Sc Reynolds 1985). The theoreti- 
cal results are then used to develop a model for the pressure-strain-rate term 
Ta — {2fp)psij. A ‘simple’ model for the pressure-strain-rate term shows that 
the ‘history effect’ (total strain) is important. 

Project C is to explore detailed processes of the pressure-strain energy transfer 
in a turbulent shear flow (Brasseur &; Lee 1987, 1988). Relationships between 
a local energy transfer event and the nearby vortical structures are studied by 
looking into instantaneous flow fields from numerical simulations of homogeneous 
shear flow (Rogers et al. 1986). For a kinematical description of the process, 
the energy transfer is classified into six classes according to the instantaneous 
directionality of the energy transfer between components. This concept of the 
energy transfer class can be applied to the statistical quantities such as the 
average value and probability density function (pdf). In this project, we are 
also interested in exploring the energy transfer across the different scales (e.g. 
spectral energy transfer) and its relationship with the intercomponent energy 
transfer. 

2. Accomplishments 

2.1. Project A: Wall Turbulence 


2.1.1. USBL Project 

The USBL project (Lee Hunt 1988) convincingly demonstrates that the 
streak generation is entirely due the high shear rate, and that the surface block- 
ing, in fact, prohibits formation of the streaks. As shown by Lee, Colonius & 
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FIGURE 2. Schematic diagram to show difference between the mechanisms 
for (a) shear-free boundary layer, where image vortex below y = 0 induces 
irrotational velocity u (B) above y = 0; and (6) uniform-shear boundary layer, 
where reduced vortex bending near the surface reduces u 2 while splat effect 
increases u 2 . 


Moin (1988), the streaks consist of eddies of large streamwise velocity fluctu- 
ations and turbulent shear stress — uv, and hence they are energy-producing 
eddies. The vorticity field associated with the streaks is found to be corrugated 
sheet vortex which has large spanwise and normal vorticity fluctuations, u> z and 
u y , but relatively weak streamwise vorticity, oj x . Near a solid surface, the cor- 
rugated vortex sheet becomes flat (less corrugated) due to the blocking effect 
of the surface, thus reducing the dominance of the streamwise fluctuations (see 
figure 2). 

Figure 3 shows the energy spectra of the streamwise velocity £’n(/C 3 ; y) as a 
function of the spanwise wavenumber and total shear f} = St at y/L = 0 
and 1 ( L is the energy integral scale at 0 = 0). While there is no evident sign 
of the streaks at the surface y = 0 (fig. 3a), the spectra away from the surface 
{y/L = l) show peaks that develop significantly with shear (fig. 3b). The 
fact that the streaks exist in sheared turbulence, independent of the presence 
of a solid boundary, but not in a shear-free boundary layer (Uzkan & Reynolds 
1967), strongly supports the assertion that the main mechanism of generating 
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FIGURE 3. Variation with shear of the spanwise one-dimensional energy 

spectra 0 ii(ks) of the streamwise fluctuation: (a) y = 0.0; (6) y = 1.0. , 

0 = 0; ,0=1; ,0 = 2] ,0 = 4; ,0-8. 
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FIGURE 4. Variation with the distance from the boundary of the mean 
streak spacing, A*. The mean streak spacing is determined from the spanwise 
wavenumber at which 0u(/c 3 ; y) peaks. 0 , 0 = 2; A , 0 = 4; X, 0 = 8. 

the streaks is the mean shear but not the wall blocking (Lee et al. 1987). 

Notice that the peaks within a narrow band of /C3 are indicative of the existence 
of the streaks in the flow. The mean spanwise spacing A, of the streaks can be 
estimated from the wavenumber at which the spectrum peaks: A* = 1 /k|. In 
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figure 4, the variation of the mean streak spacing shows an increase with y. This 
study reinforces our earlier work (Lee, Kim & Moin 1987), namely, that the 
generation mechanism of the streaks is essentially inviscid and linear. 

2.1.2. PCF Project 

The PCF and SFCF projects make use of numerical databases from direct 
simulations. We have completed numerical simulation of plane Couette flow 
(PCF) with a wall moving at a speed U w parallel to the stationary wall. The 
computation was performed on a 128 x 129 x 128 grid and a passive scalar field 
(Pr = 0.71) with constant boundary conditions was included in the computation. 
The flow Reynolds number Re based on the wall velocity U w and half the channel 
width h is 6,000, and that based on the shear velocity U T = (i/df7/dy| w ) 1//2 is 
about 200. Other major aspects of the computation are similar to those in the 
Poiseuille flow of Kim et al. (1987). The computation has been conducted for 
a period of about 100 h/U w that corresponds to about ten (10) computational 
box lengths. 

In figure 5(a), the computed mean velocity profile across the channel is com- 
pared with the experimental results conducted at different flow Reynolds num- 
bers (Re = 2, 900, Reichardt 1959; Re ~ 2 x 10 4 - 4 x 10 4 , El Telbany & Reynolds 
1980). It should be noticed that at high Reynolds numbers the velocity profile 
changes rapidly within a narrow region near the walls (|y| > 0.8), and that there 
exists a constant slope over half the channel width around the center (|y| < 0.5). 
The mean velocity gradient at the boundary dU/dy\ w grows significantly with 
the Reynolds number. The near-wall profile made dimensionless by the viscous 
units [U r and £ v = v/U T ) is in good agreement with the classical law of the wall 
(fig. 56). 

The turbulence intensities (u' + ,t/ ,+ ,u; ,+ ) scaled by the shear velocity U T 
in figure 6 also show good agreement with the experimental results at higher 
Reynolds numbers (El Telbany & Reynolds 1981). Compared with those in a 
plane Poiseuille flow at comparable Reynolds numbers (Kim et al. 1987; Kreplin 
& Eckelmann 1979), the intensities in PCF are significantly higher over most of 
the channel, except in the vicinity of the wall (y < 30) where the PCF values are 
only slightly higher. This marked contrast is a direct consequence of the con- 
stancy of total shear stress r/p = udU/dy — uv across the channel in PCF. (In 
a plane Poiseuille flow, the total shear stress has a linear profile r/|r w | = —y/h, 
and dU/dy = 0 and uv = 0 at the channel center, y = 0, by symmetry.) 

The turbulence structure exhibits strong anisotropy; even in the core region, 
u 2 : v 2 : w 2 ~ 6 : 1 : 2. The linear mean velocity profile (i.e. constant mean shear 
rate) in the core region is indicative of the existence of homogeneous turbulence, 
where the turbulence statistics are approximately constant (see figure 6). 

2.1.3. SFCF Project 

In order to achieve zero mean shear rate in the wall region of shear-free channel 
flow (SFCF), one needs to determine the speed of the moving wall U w and the 
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y/h 



Figure 5. Mean velocity profile in plane Couette flow: (a) global profile, 

U/Uit, vs. y/h; (6) near-wall profile, U + vs. y + . , present result; •, 

Reichardt (1959); H, El Telbany &; Reynolds (1980); , U + = y + ; , 

= (l/ K ) lny + + B (k = 0.4 and B = 5.5); In (6), the solid and chain-dashed 
lines are results for the lower and upper walls, respectively. 
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FIGURE 6. Near-wall profiles of the turbulence intensities in plane Couette 

flow: (u' + , t/ + ,u/ + ) vs. y + . Lines are present computation: , u ,+ ; , 

w' + ; , tv ,+ . Symbolds are from El Telbany & Reynolds (1981): O, Re 

= 19,000; 0, Re = 29,000. 

constant pressure gradient dP/dx (see figure 1), which can be obtained only by 
solving the equations. Preliminary calculation on a 64 x 129 x 64 grid has been 
carried with a developed plane Couette flow field as its initial condition; the 
value of the pressure gradient was varied at each timestep to obtain zero mean 
shear on the moving wall, and the response of the flow to the varying pressure 
gradient was found to be slow. Full-scale computation on a 128 x 129 x 128 grid 
is planned. 

2.2. Project B: Axisymmetric-strain/Dilatation RDT 

Exact expressions for turbulence structural quantities, such as the Reynolds 
stresses and vorticity correlations, were obtained in closed form by using rapid 
distortion theory (Lee 1988). Comparison with the numerical simulations (Lee 
& Reynolds 1985) at a high strain rate shows remarkable agreement for all the 
quantities considered. Differences in effects of the two opposite, axisymmetric 
strain modes (contraction and expansion) on turbulence structure were investi- 
gated in detail. For example, axisymmetric contraction produces much higher 
fluctuation in velocity and vorticity (J2jy = ttjtty and y» = WjWy, respectively) 
as well as their anisotropy (6,y = Rij/Rkk — 6»y/3 and t>,y = Vij/Vkk — <5.y/3, re- 
spectively) than does axisymmetric expansion. Figures 7(a, b) show comparison 
of anisotropy development during axisymmetric distortion. 
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FIGURE 7. Evolution of anisotropy in axisymmetric strain flows (AC, 
axisymmetric contraction; AE, axisymmetric expansion): (a) Reynolds-stress 

anisotropy ; (6) vorticity anisotropy, Lines are RDT results: AC, , 

i = 1,3 = 1; , t = 2,j = 2; AE, , * = l,j = 1; , t = 2 ,j = 2. Sym- 

bols are from numerical simulations (Lee & Reynolds 1985): AC, • , t = l,j = 1; 
O , t = 2,i = 2; AE, ■, i = 1, j =!;□,»' = 2 ,j - 2. 
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Figure 8. Test of the present model against different strain-rate cases, 
Tn/(Sql): (o) axisymmetric contraction (5 U > 0); (6) axisymmetric expansion 
< o). Lines are prediction of the present model and symbols are from 

numerical simulations (Lee & Reynolds 1985): and • , S* — 100; 

and A , S* a 10; and Q,S* a 1. 
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The effect of mean-field dilatation on turbulence is also of interest. It is shown 
that, for flows of uniform, time-dependent density p(t) at low turbulence Mach 
number, the statistical correlations, such as the ‘Reynolds stress’ Rij = u, tty , are 
related to those in the equivalent incompressible flow by a decomposition of mean 
strain-rate tensor, e.g. Rij (e) = (p/po) 2/,3j R*y( e *) where evy = f 0 Uij(t')dt' is 
the total-distortion tensor and e* } - is the corresponding incompressible distortion 
(po is the initial density). Similarly, the vorticity correlation Vii = is given 
by V{j(e) = (p/po) 4 ^ 3 Ky( e *)* ^ was a ^ so f° un( l that the rate of dilatation does 
not contribute to pressure fluctuation. 

Terms in the transport equation for were evaluated for axisymmetric strain 
flow, and existing models for the pressure-strain-rate term were examined. It 
was shown that improvements can be made for the model by incorporating struc- 
tural parameters based on the linear analysis (RDT) of axisymmetric strain 
flows. Prediction of the improved model agrees quite well with the numerical 
simulations, even in cases at lower strain rates (figure 8). This study indicates 
that a turbulence model should reflect ‘history effect, 5 since a state in nonsta- 
tionary turbulence depends not only on the instantaneous quantities (e.g. the 
Reynolds stresses) but also on the memory effect (e.g. total strain) accumulated 
during a distortion. 

2.3. Project C; Energy Transfer Process 

A detailed study of the intercomponent energy transfer processes by the 
pressure-strain-rate in homogeneous shear flow has been carried out by using a 
numerical database by Rogers et al . (1986). In the previous study (Brasseur 
& Lee 1987), it was found that t he rap id and slow parts of fluctuating pres- 
sure are uncorrelated, i.e. p7pl/(p? P?) 1 ^ 2 < 1, and their scales are widely 
separated, providing a strong justification for current modeling procedure in 
which the pressure-strain-rate term is split into the corresponding parts. It was 
shown that local intercomponent energy transfer processes are associated with 
high vorticity regions. We limit our discussion to the ‘slow 5 pressure-strain-rate 
0t/ = (2/ p) Ps s ij • 

Our recent study (Brasseur &; Lee 1988) of probability density functions 
(pdf’s) and contour plots of the pressure-strain-rate shows that the energy 
transfer processes are extremely peaky (figure 9), with high- magnitude events 
dominating low-magnitude fluctuations, as reflected by very high flatness factors 
(30-40). The concept of the energy transfer class is defined as 

C ±±=t = {0»y(x,f) |0U ^ 0,^22 ^ 0,033 ^ 0}. 

For example, C - + + denotes the class of 0<y in which energy leaves u (0n < 0) 
and enters v and w (0 22 > 0, 033 > 0). Note that this classification of the 
energy transfer is unique and there are six possible classes which are disjoint (or 
mutually exclusive). The classification has been applied to investigate details of 
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FIGURE 10. Schematic of the intercomponent energy transfer processes, 
representing contributions from all magnitudes of <f>i } - events. The schematic 
shows that the four classes C ++_ , C -++ , C and C *" are dominant. 

the direction as well as magnitude of the energy transfer processes. Examination 
of contours in an instantaneous field, pdf’s and weighted pdf’s of the pressure- 
strain-rate indicates that in the low-magnitude regions (defined as where \4> a a | < 
1.5 r.m.s. <t>aa) all six classes are important, but in the high-magnitude regions 
four classes C ++ “, C l " + , C 1 and C dominate. The contribution to 
the average slow pressure-strain-rate from the high-magnitude fluctuations is 
only 50% or less, indicating complexity of the problem. However, when summed 
over all magnitudes, the same four classes of energy transfer dominate as shown 
by the schematic in figure 10. 

3. Future Plan 


3.1. Project A: Wall Turbulence 


3.1.1. USBL Project 

We plan to look into other aspects of the blocking effect of the boundary by 
examining the two-point correlations between the vertical velocity component 
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and the streamwise and vertical components. Here, we are interested in assessing 
the degree of self-similarity of these correlations and seeking possible scaling laws 
that can be used to develop near-wall turbulence models. 

3.1.2. PCF Project 

It is anticipated to learn a great deal about how the turbulence structure 
changes across the interface between the homogeneous core region and inhomo- 
geneous wall region. We are especially interested in examining the structure of 
stationary, homogeneous, sheared turbulence in the core region. The stationar- 
ity and homogeneity in this flow are generated in a natural way, unlike in flows 
where artificial forcing techniques are used. 

3.1.3. SFCF Project 

The SFCF project complements the USBL and PCF studies in many respects. 
While the USBL and PCF projects focus on the role of the mean shear rate in 
affecting structure of turbulence in wall layer, the SFCF study is designed to 
study the kinematical nature of the presence of a boundary, i.e. the inquiry 
into how turbulence is suppressed by the boundary. We expect to learn details 
of physical mechanisms involving energy redistribution near the boundary, and 
comparison will be made between the RDT study of shear-free boundary layer, 
a special case of USBL at /? = 0 (Lee & Hunt 1988). Another impetus to carry 
out the SFCF project is that this flow allows a unique opportunity to develop 
a turbulence model for the transport terms, since, in the shear-free wall layer, 
the transport terms balance the dissipation term in the equation for turbulent 
kinetic energy. 


3.2. Project C: Energy Transfer Process 

Inspection of the energy transfer processes depicted in figure 10 poses an inter- 
esting question about the kinematic structure of the energy transfer. It appears 
that the dominant energy transfer to or from one component is provided by the 
other two components. If the off-diagonal components in 4>%j are negligible, then 
most of the energy transfer takes place within axisymmetric flow regimes. The 
answer awaits invariant analysis of the tensorial structure of the pressure-strain- 
rate (and strain rate). 

In parallel with the intercomponent energy transfer, energy transfer between 
scales in a flow with coherent vortical structures is also important. In an at- 
tempt to generate coherent structures in grid turbulence, Michard et al. (1987) 
used a grid with propellers which give rise to spectral disturbances. Their study 
showed that the initially strong anisotropy near the grid due to spectral distur- 
bances relaxes to an axisymmetric state. Numerical simulation is planned to 
investigate such a flow. We are interested in looking into how the disturbances 
are transferred between components as well as between scales. 
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Direct numerical simulation of 
compressible free shear flows 

By S. K. LELE 


Abstract 

Direct numerical simulations of compressible free shear layers in open domains 
are conducted. Compact finite-difference schemes of spectral-like accuracy are 
used for the simulations. Both temporally-growing and spatially-growing mixing 
layers are studied. The effect of intrinsic compressibility on the evolution of 
vortices is studied. The use of convective Mach number is validated. Details of 
vortex roll up and pairing are studied. A simple explanation of the stabilizing 
effect of compressibility is offered. Acoustic radiation from vortex roll up, pairing 
and shape oscillations is studied and quantified. 

1. Introduction 

Direct numerical simulations have been used to study several incompressible 
turbulent flows. These studies have ranged from the simulations of homoge- 
neous turbulence (Rogallo, 1981, Lee and Reynolds, 1985, and Rogers, Moin and 
Reynolds, 1986), to the turbulence in channel flows, (Moser and Moin, 1987, and 
Kim, Moin and Moser, 1987), boundary layers, (Spalart, 1988), and free shear 
flows, (Corcos and Sherman, 1984, Corcos and Lin, 1984, Riley and Metcalfe, 
1980, Metcalfe et. al., 1987, Lowery and Reynolds, 1986, and Sandham and 
Reynolds, 1987). The results obtained from the simulations have been exten- 
sively compared with the experimental measurements (Kim, Moin and Moser, 
1987, Spalart, 1988) and are being used as databases for developing and testing 
turbulence models (Mansour, Kim and Moin, 1988, and Center for Turbulence 
Research, 1987), studying the structure and dynamics of organized motions in 
turbulent flows (Center for Turbulence Research, 1987 and 1988), and test- 
ing the accuracy of experimental techniques by Moin and Spalart (1987). On 
the other hand, the study of compressible turbulent flows is relatively less ma- 
ture. Most simulations of compressible flows have been limited to the simula- 
tion of the large scales of the flow (Winkler et. al., 1987, Smarr et. al., 1984, 
Woodward, 1984, and Boris, 1988), with some notable exceptions (Feiereisen, 
Reynolds and Ferziger, 1981, and Passot and Pouquet, 1987). Comparisons with 
experimental measurements are less detailed in part due to the ‘large scale’ and 
two-dimensional nature of most simulations and lack of experimental data on 
turbulent fluctuations in compressible flows. The renewed interest in the high 
speed flows, as well as the need to develop turbulence models which incorporate 
the compressibility effects provide the motivation for the present work. This 
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paper provides a brief account of results obtained from direct simulations of 
compressible free shear layers. The simulations resolve all the relevant scales of 
motion and in particular do not employ any form of filtering or artificial dissipa- 
tion. Simulations are conducted for both two-dimensional and three-dimensional 
flows. Most results presented are based on the two-dimensional simulations. 


2. Numerical Method 

The numerical method used in the study is a compact finite difference scheme 
of spectral-like accuracy. The details of the numerical scheme have been de- 
scribed in detail elsewhere (Lele, 1988a). The full unsteady equations for mass, 
momentum and energy conservation are solved for an ideal gas. The Newtonian 
viscosity law (zero bulk viscosity), and Fourier law of heat conduction are used. 
Variation of viscosity and thermal conductivity with temperature is permit- 
ted. Conservation equation for of a passive scalar (mass fraction) is also solved. 
Time advancement is carried out by a compact storage third order Runge-Kutta 
scheme developed by Wray (1987). 

The main features of the finite-difference scheme are summarized. The fi- 
nite difference approximation to the first derivative ^ (z*) at the node * is 
evaluated by solving a tridiagonal system of the form: 


ifi+2 — fi - 2 , fi+1 ~ fi-l 
a/i-i + fi + af i+1 = b - + a- 


4h 


2h 


( 1 . 1 ) 


The relations between the coefficients a,b and a are derived by matching the 
Taylor series coefficients of various order. With 


°= |(a + 2 ),6 = |(4a- i) (1.2) 

a family of fourth order schemes is obtained. It may be noted that as a — ► 0 
this family merges into the well known fourth order central difference scheme. 
Similarly for a = 4 the classical Pade’ scheme is recovered. Furthermore for 
a = ^ the leading order truncation error coefficient vanishes and the scheme 
is formally sixth order accurate. Most of the simulations presented here use 
this sixth order scheme. The spectral-like accuracy of the scheme follows from 
the nature of the dispersive errors associated with (l.l). It may be shown that 
compared to the traditional finite difference schemes the scheme (l.l) reduces 
the dispersive errors over a wider band of the length scales represented on the 
grid. This characterization of the dispersive errors is presented in figure-1. The 
straight line on the figure represents a spectrally accurate scheme (i.e. exact for 
all the wavenumbers represented on the mesh). The improved representation of 
the shorter spatial scales by the present scheme is evident from the figure. 

If the dependent variables are periodic then the system of relations (l.l) writ- 
ten for each node can be solved together as a linear system of equations. The 
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general non-periodic case requires additional relations appropriate for the near 
boundary nodes. These are described in Lele (1988a). The relation (1.1) along 
with a mathematically defined mapping between a non-uniform physical mesh 
and a uniform computational mesh provides derivatives on a non-uniform mesh. 
The second derivatives are evaluated by solving a system similar to (l.l), viz. 

a fi-l + fi + “/» + 1 = 

l /»+2 — 2/i + /„•— 2 , fi+l ~ 2fi + fi-1 
b 4 *2 + “ • (»•») 

where /" represents the finite difference approximation to the second derivative 
at node i. With 

o=^(l-a),6=i(-l + 10 a), (1.4) 

a one parameter family of fourth order schemes is obtained. Again as a — ► 0 
this family coincides with the well known fourth order central difference scheme. 
For a = the classical Pade’ scheme is recovered. For a = a sixth order 
tridiagonal scheme is obtained. This scheme with 




3 

II 


(1.5) 


is used for the simulations presented here. It may also be noted that the schemes 
(1.3) provide an accurate evaluation of the second derivative over a wide range 
of length scales. The error associated with the second derivative evaluation for 
a variety of schemes is shown on figure-2. The spectrally accurate evaluation 
corresponds to the parabola on the figure. The improvement of the present 
scheme in representing the shorter scales is again evident from this figure. 

Simulations are conducted for both spatially-evolving and temporally-evolving 
mixing layers. In both cases non-reflecting boundary conditions (Thompson, 
1987) are employed at the top and bottom boundaries of the computational do- 
main. In the spatially-evolving case non-reflecting boundary conditions are also 
used at the inflow and outflow boundaries. The computational grid is uniformly 
spaced in x (mean flow direction) and z (spanwise direction) and non-uniform 
in y (transverse to the mean flow). Typically a hyperbolic tangent mapping is 
used for the non-uniform mesh (in y) with the maximum grid spacing about 2-4 
times the minimum grid spacing. 


3. Mixing Layer Flows 

Direct numerical simulations of mixing layer forming between two streams 
flowing with speeds Ui (faster stream) and were conducted. The two streams 
were matched in static pressure but, in general, had different fluid densities. Such 
two stream mixing layers were studied by Brown and Roshko (1974) and subse- 
quently by other investigators. They noted that while the density ratio of the 
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two streams had a noticeable effect on the spreading rate of the mixing layer it 
was not enough to explain the slow spreading rate of supersonic mixing layers. 
Papamoschou and Roshko (1986) (referred to as PR from hereon) showed that 
this compressibility effect on the spreading rate could be parameterized in terms 
of a single Mach number, the Convective Mach Number M c , viz. the Mach num- 
ber defined in a frame of reference moving with the dominant eddies of the flow. 
They also derived an expression for M c in terms of the flow conditions of the 
two streams. Bogdanoff (1983) also suggested the same ‘intrinsic Mach number’ 
by a different argument. Recent experiments (Papamoschou, 1988, and Samimy 
and Elliott, 1988) stability analyses (Sandham and Reynolds, 1989, Ragab and 
Wu, 1988, Zhuang et. al., 1988) and numerical simulations (Soetrisno et. al., 
1988, Lele 1988b) have provided further support to the this notion. Numeri- 
cal simulations are used to test the notion of the convective Mach number in a 
fundamental way. Simulations were conducted for three types of mixing layer 
flows: A) Mixing of streams of equal entropy, B) Mixing of streams of equal 
stagnation enthalpy, and C) Mixing of streams of equal Mach number. These 
different classes of mixing layer flows allowed us to vary the Mach number of 
the individual streams while keeping the convective Mach number fixed. It is 
thus possible to look for the influence of factors other than the convective Mach 
number. The spatially evolving simulations presented here were forced at the 
inflow by adding a v-velocity disturbance velocity of small amplitude (typically 
with a peak of l%)to the inflow profiles. This disturbance was confined to the 
shear layer by means of a Gaussian shape function. It is possible to refine the 
forcing by using the linear stability eigenfunctions of the inflow profiles. Such 
a approach was used by Sandham and Reynolds (1987) in simulating incom- 
pressible mixing layers. The inflow disturbance contained the most unstable 
frequency and its first two subharmonics. 

In the temporally evolving simulations the initial condition provided the source 
of disturbance. Small amplitude incompressible disturbances were added to the 
tangent hyperbolic mean flow profile. Two groups of simulations were con- 
ducted. l) Organized initial conditions - the disturbance contained a small 
number (typically two) of fourier modes with prescribed phase relations, and 
2) Random initial conditions - the disturbance had a white noise wavenumber 
spectrum with very small (typically 10 -4 or smaller) amplitude. Both 2-d and 
3-d simulations were conducted with the random initial conditions. 

The simulations had specific heat ratio, 7 = 1.4, Prandtl number, Pr = 3/4, 
and Reynolds number R = P ^ - 1 in the range 100-500, where subscripts 

1 and 2 refer to the high speed and low speed streams, respectively, and 6 U , 0 is 
the initial vorticity thickness. 

4. Intrinsic compressibility 

The definition of a convective Mach number requires one to specify the speed of 
propagation of the dominant eddies. This speed of propagation was determined 
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in the numerical simulations from plots of the location of local pressure ex- 
tremas against time. The pressure maxima correspond to the stagnation points 
between the vortices, and the pressure minima correspond to the vortex centers. 
The propagation speed is taken as the slope of these trajectories. An example 
of such a plot is shown in figure 3 for a class-A flow. For this particular case the 
Mach numbers of the two streams are 2.0 and 1.2 and the inflow disturbance 
contains the fundamental and its first two subharmonics. The trajectories of the 
stagnation points are shown in figure 3a. The location of the vortex centers are 
displayed in figure 3b. Two generations of pairing events can be seen in these 
figures. It may be noted that even with the pairing events the stagnation points 
move with a relatively uniform speed. During the pairing events the downstream 
vortex slows down while the upstream vortex speeds up. As the vortices pair 
the stagnation point in their middle is lost. It is also seen that the lower subhar- 
monics tend to modulate the location of upstream pairing events. Varying the 
frequency of the fundamental disturbance did not alter the propagation speed 
of the vortices. 

PR presented a formula for the vortex propagation speed based on the phys- 
ical argument that the static pressure at the stagnation point is related to the 
stagnation pressures in the two streams by the isentropic relation. This convec- 
tion speed formula was tested against the numerical simulation results at three 
different velocity ratios and several convective Mach numbers. The agreement 
between the observed speed and the formula is excellent in all the simulations 
(discrepancy of 1 - 2 %). It may be noted that when the specific heat ratio of 
the two gases is the same, the PR formula reduces to the formula proposed by 
Dimotakis (1984) for the propagation speed of vortices in incompressible shear 
layers. Thus fluid compressibility appears to have very little influence on the 
propagation of two-dimensional vortices. 

The spatial growth rate of the mixing layer depends both on the growth mea- 
sured in the convected frame of reference (temporal growth) and the convection 
speed of the vortices. Having verified that the convection speed is indepen- 
dent of the compressibility, compressibility effects in the spatial growth of the 
layer are studied. Figure 4 shows an example of the downstream evolution of 
three measures of mixing layer thickness. These are: 1) Thickness based on the 
Reynolds-averaged mean flow profile, 2) Thickness based on the Favre-averaged 
mean flow profile, and 3) Thickness based on the mean potential vorticity profile. 
For the example displayed, the layer was forced at the most unstable frequency 
with a 5 % forcing amplitude. The three cases (one from each class) have M c 
of about 0.4 and the same velocity ratio ^ of 0.6. The thickness is normalized 
by the initial thickness and the downstream distance is normalized by 10 6 Uo . It 
may be seen when M c is the same the spatial growth rate for the three cases is 
similar but not identical. The differences that remain arise from the dependence 
of the nondimensional convection speed ^ on the density ratio ^ . Case C has 
the slowest propagation speed giving the fastest spatial growth, while case B has 
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the fastest propagation speed and the slowest growth. A visual picture of the 
spatially evolving mixing layers is presented in figure-5. Contours of vorticity are 
plotted for three simulations with M c about 0.4. The simulations were forced 
with the fundamental and two subharmonics. Vortex pairing events noted in 
figure-3 can be visually seen in these contours. It may also be noted that the 
vortices are most evolved in case C (slowest propagation speed). 

Results obtained from the vortex roll up in temporally growing mixing layers 
are presented to further verify this density ratio effect on The temporal 

simulations are conducted in a frame of reference moving with U c = , 

ax and 02 being the sound speeds in the two streams. In this way the density 
ratio effect is removed. The initial amplitude in all the temporal simulations is 1 
%. In figure 6 two examples are shown with M c of about 0.4. The first case has 
streams with equal static temperature, i.e. T 2 = Tx (class A) and the second case 
has streams with equal total temperature, i.e. h 01 = ho 2 (class B). Thickness 
measures (l) and (2) defined above are displayed against time normalized by 
. It may be noted that the time histories are almost identical. 

Experiments have documented the reduction in the growth rate of compress- 
ible shear layers as M c is increased (PR). This reduction is most easily under- 
stood in terms of the behavior of the temporally-evolving layer. In figure 7 we 
show the time history of thickness measure (2) for several different M c . The 
stabilizing influence of M c is clearly seen. It has been suggested that the slow 
growth arises primarily due to the reduced linear instability growth rate (PR 
and Sandham and Reynolds, 1989). 

In the next section we present a physical argument for this stabilizing effect. 


5. Evolution of the vorticity field 

The linear instability process leads to a reorganization of the vorticity field. 
The inviscid vorticity equation for the two dimensional flow may be written as: 


du ,, dw , __ . du du _ 


_ Vox Vp 
g+ - „ - - 


showing that the rate of change of vorticity observed by an observer moving 
at speed U c is due to three effects. 1) Advection relative to the observer (first 
two terms on the r.h.s.), 2) Change due to compression of fluid elements (third 
term), and 3) Baroclinic change (fourth term). The individual terms for a class 
A mixing layer at M c of 0.4 are shown in figure 8. It may be noted that the 
most important term causing the vorticity redistribution is the advection term. 
Advection may be seen as moving vorticity away from the stagnation region and 
bringing it to the vortex centers. This is precisely the incompressible instability 
mechanism causing the shear layer instability. The compression term, while not 
dominant at M c = 0.4, provides clues on how the stabilizing effect associated 
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with compressibility arises. A simple explanation of the stabilizing effect of com- 
pressibility is presented below. In figure 10 a schematic picture of the flow field 
around the vortices is presented in a frame of reference moving with speed U c , 
the speed of the vortices. As shown earlier the ‘stagnation points’ between the 
vortices also move at the same speed and thus are stationary in this frame. Fluid 
moving away from the stagnation points expands as it accelerates towards the 
vortex center. Past this point the flow compresses and decelerates towards the 
other stagnation point. This expansion and compression process may be seen 
in the characteristic quadrupole pattern of the dilatation field. The expansion 
and compression cycle has the effect of increasing the vorticity near the stagna- 
tion points and reducing it near the vortex center. This is exactly opposite of 
the redistribution arising from the advection term (the driving cause of insta- 
bility). The numerical simulations show that as M c is increased the vorticity 
compression effect becomes comparable in magnitude to the advection term. 

In class A flows, the baroclinic change is an order of magnitude smaller than 
the other terms. In class B and C flows, the density of the two streams are 
unequal and this generates strong baroclinic torques. This is essentially an 
incompressible effect and arises whenever a density gradient exists across the 
shear layer. As the roll up proceeds this density interface remains sharp at the 
‘braids’ and the pressure maximum (at the stagnation point in the convected 
frame) produces regions with dynamically significant baroclinic torque. For 
class B flows, the density of the low speed stream is always lower, thus the 
baroclinic torques tend to enhance the vorticity in the lower part of the braid 
and cause suppression of vorticity in the upper part of the braid. For class C 
flows, the density of the high speed stream is lower and a reversal of the situation 
just described takes place. The vortices evolving in class B and class C flows 
have a layered vorticity distribution. The baroclinic effect strongly modifies the 
vorticity distribution in the vortices. This mechanism of generation of baroclinic 
torques is shown in figure 9 for a class C flow. The net circulation of the vortices, 
however, remains relatively uninfluenced by this redistribution. 

6. Flows with eddy shocklets 

For M c around 0.7 and larger, the flow fields develop eddy shocklets. These 
shocks remain attached to the vortices and travel with them. For M c = 0.7 
the shocks arise during the vortex pairing events, but the flow field is other- 
wise shock free. For M c — 0.8 the shocklets are produced during the roll up. 
An example of flow fields with eddy shocklets is displayed in figure 11. Just 
as observed at lower M c the flow accelerates and decelerates around the vor- 
tices. Now, however, there are local regions where the relative Mach number 
of unity is exceeded, and near the vortex this supersonic flow slows down by 
first going through a shock, becoming subsonic, and then decelerating further 
by compressing towards the stagnation points. The regions of expansion are 
now more spread out and the regions of compression more compact (essentially 
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within the eddy shocklets). The increase of the vorticity due to the compression 
in the shock, as well as the increase in entropy is clearly observed. This flow 
field is remarkably similar to typical transonic flow past an airfoil. The induced 
velocity pattern associated with the fluid expansion and compression opposes 
the entrainment velocity induced by the clumped vorticity field. This may cause 
a further reduction in the entrainment of fluid into the layer. 

7. Acoustic radiation from vortex evolution 

Numerical simulations have been used to study the acoustic radiation arising 
from the unsteady processes of vortex roll up, pairing, shredding, and shape os- 
cillations. The acoustic efficiency for these processes has been studied at several 
Mach numbers. At low Mach number less than 0.1% of the energy extracted 
from the mean flow is radiated. This fraction increases to 1-2 % at M c of 0.6. 
A full account of these studies will be presented elsewhere (Lele and Ho, 1989). 
Here we present one example to illustrate the behavior. The example chosen 
corresponds to a temporally-growing mixing layer. The computation uses non- 
reflecting boundary conditions based on characteristic variables. These allow 
the acoustic waves to leave the computational domain. The example discussed 
here has M c = 0.6 , a value for which we find significant departures from the 
low Mach number aero-acoustic theory. The computational domain contains 
two wavelengths of the most unstable disturbance. The layer first rolls up to 
form two vortices which later pair to form a larger vortex. After the pairing the 
vortex continues to undergo shape oscillations for several cycles. This behavior 
is illustrated in figure 12 where the time history of layer thickness is plotted. 
The corresponding time history of pressure, velocity component normal to the 
layer, and the acoustic flux of energy leaving the domain near the bottom are 
also shown in this figure. Pressure is normalized as velocity as - } ^r, 

and acoustic flux in terms of these normalized values. The roll up, pairing and 
nutation processes can be identified in these figures. The roll up and pairing 
cause the layer thickness to increase and generate a compression wave. Subse- 
quent to pairing, the nutation of the vortices causes a periodic energy exchange 
between the vortices and the mean flow. This process generates a series of com- 
pression and expansion waves which carry acoustic energy away from the layer. 
The numerical simulations show that the nutation frequency is close to ^ , where 
u> is the peak vorticity. 

In figure-13 the acoustic radiation for three different cases are compared. All 
case were for M c = 0.4, but differed in the initial disturbance. Case-1 con- 
tained only the fundamental disturbance (case-1), case-2 only the subharmonic 
disturbance, and case-3 contained both the fundamental and the subharmonic. 
It is evident (as anticipated by Laufer (1974)) that in case-3 the vortex merger 
produces the strongest acoustic radiation (8 times that in case-1 and twice that 
of case-2). In the time history of the acoustic signal the signals arising from the 
rollup and pairing can be separately identified. The acoustic signal from rollup 
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contains a compression wave followed by a more spread out expansion wave 
while the acoustic signal from pairing generates more symmetric compression 
and expansion waves. 

8. Simulations with random initial conditions 

Temporally evolving mixing layer flows were simulated starting from very low 
level (peak amplitude 10~ 4 ) white noise (equal amplitude to all fourier modes 
with uncorrelated phases). Both two-dimensional and three-dimensional simu- 
lations were conducted. The computational box was chosen to be 12 times the 
most unstable wavelength (based on the initial profile). The initial evolution of 
the disturbances corresponded to the predictions of the linear stability analysis. 
In figure-14 the exponential growth rates obtained from two such simulations are 
shown. The two cases correspond to M c = 0.4 and M c = 0.8. The wavenumber 
and the growth rate are scaled with the local values of the vorticity thickness. 
The time scale is normalized by u*-u o in this figure. During this phase of ex- 
ponential amplification the layer grew in a laminar fashion. As the disturbances 
became nonlinear the harmonics of the linearly unstable disturbances began to 
amplify. With significant vorticity clumping the growth rates became larger than 
the ‘laminar instability’ estimates. The nonlinear process generated a broadband 
spectrum. In this regime a linear spreading of the mixing layer was observed. 
In figure-15 data from three simulations are presented. It may be seen that for 
M c = 0.4 and M c = 0.6 a constant spreading rate was observed, with slower 
growth for M c = 0.6. It was verified that for M c lower than 0.4 the growth rate 
was not further increased. The temporal growth rate obtained for M c = 0.4 is 
close to 0.1 a value consistent with previous incompressible simulations Lesieur 
(1987). For M c = 0.8 the growth was slow and dominated by a laminar process. 
Later in time (outside the time range shown) vortices formed in the layer but 
the computational domain contained too few of them to constitute an adequate 
statistical sample. The convective Mach number concept was also tested with 
random initial conditions. In figure-16 the time evolution of two run with the 
same M c of 0.4 but different density ratios are compared. The temporal growth 
rates were quite close (though not identical). 

In figure-17 snapshots of the vorticity field are displayed for the case with 
M c = 0.4. The frames are equally separated in time. It may be seen that 
the linear instability process selects a length scale for the vortices. This length 
scale does not correspond to the most unstable wavelength of the initial profile. 
This is because the initial disturbance level was very low, allowing the layer to 
viscously thicken [Res Mo = 400). The length scale chosen by the flow corresponds 
to the most unstable wavelength at a later time (end of the laminar spreading 
regime in fig. 15). After the vortices form they undergo the merging instability. 
Since the vortices had a phase jitter the merging is also randomized. Figure-18 
presents similar snapshots for the mixing layer at M c = 0.6. It maybe noted 
that the length scale selected by the flow is larger. This is consistent with the 
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linear instability behavior. The phase jitter in the vortices is again visible. Such 
phase jitter allows the turbulence statistics to evolve in a self-similar manner. 
Figure- 19 presents some of these statistics at four instants in the evolution (from 
the linearly spreading regime). The y-coordinate in the plots is scaled with the 
instantaneous vorticity thickness and the fluctuations are scaled by half the 
velocity difference across the layer. The mean velocity profiles are closely self- 
similar. The Reynolds stress is also close to being self-similar. The r.m.s. u and 
r.m.s. v (titled as u prime and v prime) have larger departure from self-similarity. 
Later in the evolution as the second pairing takes place the computational box 
contains too few vortices (‘independent’ samples) and the layer departs from 
self-similarity. As noted earlier the nonlinear process generate higher harmonics 
and create a broadband spectra. In figure-20 spectra of the fluctuating fields 
are plotted at various times. The spectra were obtained by fourier decomposing 
the fields in the periodic direction. The wavenumber was normalized with the 
initial vorticity thickness. The plots show the spectra obtained by integrating 
the energy associated with each fourier mode across the layer. The simulation 
start from a low level white noise. Linear instability selectively amplifies the 
unstable modes. Later in time the spectra fill up in the harmonics and finally 
become broadband. It is possible to see a shift in the peak of the spectra to low 
wavenumbers (since the layer grows in time) . Even at late times the energy level 
of the broadband fluctuations is quite low. Previous studies of incompressible 
mixing layer (Lesieur, 1987) have suggested a -4 slope in the broad band portion. 
Such a spectral slope may be seen in the present simulations. Interestingly the 
low wavenumber portion of the pressure spectrum is very much like the velocity 
spectrum, but at higher wavenumber the pressure spectrum drops more rapidly. 

This work was presented at the AIAA meeting in Reno, 1989. It is available 
as AIAA paper AIAA-89-0374 from AIAA. 

REFERENCES 

BOGDANOFF, D. W. 1983 AIAA J., 21, 926-927. 

BORIS, J. P. 1988 in Proceedings of the conference on ‘Physics of Compressible 
Turbulent Mixing’, Princeton, (see also Oran, E.) 

Brown, G. L. AND ROSHKO, A. 1974 J. Fluid Mech., 64, 775-816. 

CENTER FOR Turbulence Research 1987 Studying Turbulence using Nu- 
merical Simulation Databases, Proceedings of the 1987 Summer Program, 
Report CTR-S87, Stamford University. 

CENTER FOR TURBULENCE RESEARCH 1988 Studying Turbulence using Nu- 
merical Simulation Databases - II, Proceedings of the 1988 Summer Pro- 
gram, Report CTR-S88, Stanford University. 

CORCOS, G. M. AND SHERMAN, F. S. 1984 J. Fluid Mech., 139, 29-65. 



Direct numerical simulation of compressible free shear flows 


89 


CORCOS, G. M. AND LIN, S. J. 1984 J. Fluid Mech., 139, 67-95. 

DIMOTAKIS, P. E. 1984 AIAA paper, AIAA-84-0368. 

Feiereisen, W. j., Reynolds, W. C. and Ferziger, J. H. 1981 
Dept, of Mech. Engrg., Stanford University, Report No. TF-13. 

KIM, J., MOIN, P. AND Moser, R. D. 1987 J. Fluid Mech., 177, 133-166. 

LAUFER, J. 1974 Omaggio a Carlo Ferrari , 451-464. 

LEE, M. J. AND REYNOLDS, W. C. 1985 Dept, of Mech. Engrg., Stanford 
University, Report No. TF-24. 

LELE, S. K. 1988a submitted to J. Comput. Phys. 

LELE, S. K. 1988b in Proceedings of the conference on ‘Physics of Compress- 
ible Turbulent Mixing’, Princeton. 

LELE, S. K. AND HO, C. M. 1989 in preparation. 

LESIEUR, M. 1987 Turbulence in fluids, M. Nijhoff Publishers, Chapter 8, p253. 

LOWERY, P. S. AND REYNOLDS, W. C. 1986 Dept, of Mech. Engrg., Stan- 
ford University, Report No. TF-26. 

MANSOUR N. N., Kim, J. AND Moin, P. 1988 J. Fluid Mech., 194, 15-44. 

METCALFE ET. AL. 1987 J. Fluid Mech., 184, 207-243. 

Moin, P. AND SPALART, P. R. 1987 NASA TM-100022. 

MOSER, R. D. AND MOIN, P. 1987 J. Fluid Mech., 175,479-510. 

PAPAMOSCHOU, D. AND ROSHKO, A. 1986 AIAA paper, AIAA-86-0162. 

PAPAMOSCHOU, D. 1988 in Proceedings of the conference on ‘Physics of Com- 
pressible Turbulent Mixing’, Princeton. 

PASSOT, T. AND POUQUET, A. 1987 J. Fluid Mech., 181, 441-466. 

RAGAB, S. A. AND Wu, J. L. 1988 AIAA paper, AIAA-88-0038. 

Riley, J. J. and Metcalfe, R. W. 1980 AIAA paper, 

ROGALLO, R. 1981 NASA TM-81315. 

Rogers, M. M., Moin, P. and Reynolds, W. C. 1986 Dept, of Mech. 
Engrg., Stanford University, Report No. TF-25. 

SAMIMY, M. AND Elliott, G. S. 1988 AIAA paper 88-3054A. 

SANDHAM, N. D. AND Reynolds, W. C. 1989 AIAA paper, AIAA-89- 
0371. 

SANDHAM, N. D. AND REYNOLDS, W. C. 1987 in Turbulent Shear Flows - 
VI , Toulouse, France. 

Smarr, L. l., Norman, M. l. and Winkler, K-H. A. 1984 Physica 
D, 12, 83-106. 



90 


S. K. Lett 


SOETRISNO, M., EBERHARDT, S., RILEY, J. J., AND MCMURTY, P. 
1988 AIAA paper, AIAA-88-3676. 

SPALART, P. R. 1988 J. Fluid Mech ., 187, 61-98. 

THOMPSON, K. W. 1987 J. Comput. Phys., 68, 1-24. 

WINKLER ET. AL. 1987 Phys. Today , October, 28-37. 

WOODWARD, P. 1984 L. L. L. preprint UCRL-90662. 

WRAY, A. A. 1987 submitted to J. Comput. Phys. 

ZHUANG, M., KUBOTA, T. AND DIMOTAKIS P. 1988 in Proceedings of the 
l at National Fluid Dynamics Conference, Cincinnati. 88-3583-CP. 



Modified Wavenumber Modified Wavenumber 

0.0 1.0 2.0 3 0 1.0 5.0 6.0 7.0 B.O 9.0 0.0 0.S 1.0 t.S 2.0 


Direct numerical simulation of compressible free shear flops 


91 



Figure 1 Modified wavenumber analysis for the ap- 
proximations to the first derivative. 





Figure 3 Trajectories of pressure extrema in a spa- 
tially evolving mixing layer. Pressure max- 
ima shown in a) correspond to the stagna- 
tion points in the convected frame and pres- 
sure minima shown in b) correspond to the 
vortex centers. 


Figure 2 Modified wavenumber analysis for the ap- 
proximations to the second derivative. 
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Figure 4 Stream wise evolution of mixing layer thick- 
ness measures. Ail cases ha-e ft of 0.6, 
Me about 0.4 and the same inflow forcing. 
The difference in the growth rates is due to 
the different density ratios. A) J\ = Tj, B) 
To, = T 0 „ C) Mx = 




Figure 6 Temporal evolution of mixing layer thickness 
measures. Both cases hare the same M e = 
0.38, A) Ti = Tj, B) To, = To,. 



Figure 7 Temporal evolution of mixing layer thickness 
measures for different M c . All cases are from 
class A : a) Af e = 0.2, b) M t = 0.38, c) 
M e = 0.6, d) = 0.8. 
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Figure 8 Vorticity evolution for class A mixing 
layer. M\ = 2.0, M 2 - 1.2, Af c = 
0.4. The panels are : vorticity field, 
advection term in moving frame (- 
0.91 to 0.76), vorticity compression 
term (-0.22 to 0.25), and baroclinic 
term (-0.03 to 0.04), respectively. 
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Figure 9 Vorticity evolution for class C mix- 
ing layer. M\ — 1.5, M 2 = 1.5, 
M c = 0.38. The baroclinic terms 
are comparable to the vorticity ad- 
vection and produce the layered vor- 
ticity profiles. The panels are : vor- 
ticity field, pressure field, density field, 
and the baroclinic term, respectively. 

E : Expansion 
<7. U > 0 

C ; Compression 

V U< 0 


Figure 10 A schematic of the flow field in a frame of reference moving with the vortices. 










Figure 11 Vortex evolution with eddy shocklets. The 
example shown is for class A mixing layer. 

= 4 . 0 , M 7 = 2 . 4 , M e - 0 . 8 . The pan- 
da present the /orticity field, pressure field, 
density field, and temperature Held, respec- 
tively. 
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Figure 12 Acoustic radiation from vortex pairing and 
shape oscillation. M € = 0.6, Time history of 
mixing layer thickness measures, normalised 
pressure, normalised v velocity, and normal- 
ised acoustic flux at the bottom boundary. 



Figure 13 Acoustic radiation from vortex evolution. All 
cases have M t = 0.4. Case-1 (solid) has 
only the fundamental disturbance, 2) only 
the subharmonic disturbance (dashed) and 
3) both togather (chain). The normalised 
pressure signal in the far-fteld and normal- 
ised acoustic flux at the top boundary are 
shown. 
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Figwe 15 Mixing layer growth for different M e . The 
curves correspond to M € = 0.4, 0.5 and 0.8. 



Time 

Figure 14 Growth rates for different fourier modes for 
M c = 0.4 (solid) and Ai c = 0.8 (dashed). 
During this phase laminar spreading is ob- 
served. 



Figure 1« Mixing layer growth at M € = 0.4. A) 7\ = 
Tj. B) T 0l = T 0l . 
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Direct simulation of compressible reacting flows 

By T. J. POINSOT 


Summary 

A research program for direct numerical simulations of compressible reacting 
flows is described. Two main research subjects are proposed: the effect of pres- 
sure waves on turbulent combustion and the use of direct simulation methods 
to validate flamelet models for turbulent combustion. The interest of a com- 
pressible code to study turbulent combustion is emphasized through examples 
of reacting shear layer and combustion instabilities studies. The choice of exper- 
imental data to compare with direct simulation results is discussed. A tentative 
program is given and the computation cases to use are described as well as the 
code validation runs. 

1. Introduction 

A considerable part of progress to be made in future years in the field of turbu- 
lent combustion studies might be achieved through direct simulation methods. 
These methods have already demonstrated their possibilities in the case of non- 
reacting flows, and there is little doubt that they will be as powerful in the case 
of reacting flows. However, such progress will require considerable improvements 
of existing direct simulation codes to account for specific phenomena occurring 
in reacting flows. One of the most important issues in this field is to take into 
account the compressibility effects. This should be done at different levels by 
incorporating the following effects: 

1- density variations generated by heat release 

2- acoustic waves 

S- strong compressibility effects due to high relative Mach numbers 

These three steps represent growing complexity. It is rather obvious that den- 
sity variations generated by heat release strongly modify the velocity field and 
the strain rates, changing the turbulent field and, therefore, the combustion rate 
itself. Therefore, direct simulations aimed at studying cold flames (essentially 
flames with constant density) can not predict the behavior of real flames. This 
does not mean that these simulations are not interesting: they retain many fea- 
tures of real flames like chemical and diffusive mechanisms. Building a code 
for direct simulation of flames with non constant density does not necessarily 
mean using a compressible code. Random vortex methods as well as spectral 
methods can provide such capacity without taking compressibility effects into 
account (Ghoniem and Sethian 1987, Riley et al 1986, Givi and Metcalfe 1986). 
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Nevertheless, ability to include compressibility is highly desirable because it ful- 
fills not only criterion 1 but also the two others: computing acoustic waves and 
mean pressure gradients. 

Experimental results suggest that acoustic waves are of primary importance 
for combustion problems (Poinsot et al 1986, Darabiha et al 1985) but also for 
non-reacting flows. Non-reacting shear layers exhibit phenomena where acous- 
tic waves play an important role by introducing a feedback loop through which 
phenomena occurring downstream (for example a vortex coalescence) can send 
information upstream (to the vortex shedding region) and create a locking phe- 
nomenon (Ho and Huerre 1984, Ho and Huang 1982). This locking may be weak 
or strong. We will call it weak when it does not lead to an instability where the 
acoustic field becomes intense. This is the case in most ducted shear layers. It 
should not be concluded that weak coupling means no effect of acoustic waves : 
even if the acoustic field is not intense, vortex shedding can be very sensitive to 
pressure waves, and this point is seldom taken into account in direct simulation 
methods. 

In non-reacting shear layers the coupling may become strong as, for example, 
in the case of the edgetone experiment. Placing an obstacle in the shear layer can 
amplify the acoustic feedback and lead to the production of a highly coherent 
and intense sound generated by the coupling between coherent structures in the 
shear layer, their impact on the obstacle, and the acoustic waves in the flow (Ho 
and Nosseir 1981, Tang and Rockwell 1983, Knisely and Rockwell 1982). 

When combustion takes place, the coupling is even enhanced by the reac- 
tion taking place in the flow structures. Experimentally, this coupling is strong 
enough in many situations to induce very large oscillations of all the flow pa- 
rameters and lead to complete extinction by blow off or flash back of the flames. 
Although it does not always lead to such extreme instabilities, we have to suppose 
that this coupling can be an important part of the global behavior of turbulent 
reacting flows and, therefore, incorporate acoustic waves in our model. 

It is possible to ask whether we really need direct simulation codes to study 
interactions between combustion and acoustic waves and whether we could sim- 
plify the problem and for example assume that acoustic waves and turbulence 
act independently on the flow. Acoustic waves in usual systems have large wave- 
lengths compared to the reaction zone thickness, and we might try to disconnect 
the effects of pressure waves and those of turbulence. Such attempts have been 
made to compute combustion instabilities in the ‘Thin Flame* models proposed 
by Yang and Culick (1986) or Poinsot and Candel (1988). In these models the 
interaction between combustion and turbulence is represented simply by a tur- 
bulent flame speed while all interactions between the flame movements and the 
acoustic field are explicitly computed. Poinsot and Candel (1988) show that 
these models are unable to predict any turbulent combustion instability if the 
turbulent flame speed does not depend on the local turbulence or on the local 
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pressure. This shows that pressure waves and turbulence have strong interac- 
tions and that this coupling must be explicitly taken into account. This requires 
the use of a direct simulation code (for turbulence estimation) which must be 
compressible (for acoustic waves computation). 

The last case requiring compressible computations is clearly obtained for flows 
exhibiting high relative Mach numbers such as supersonic reacting flows. Com- 
pressibility is not limited to the flow perturbations (as it was for acoustic waves) 
in these cases but is a basic characteristic of the mean flow itself. Compressible 
calculations are obviously required for this situation. 

Looking at all these cases is far beyond the scope of this work. Therefore, 
I have chosen to concentrate this study on two generic cases of compressible 
reacting flows which offer a very wide range of interests: 

> Effect of acoustic waves on reacting and non-reacting shear flows 

> Validation of flamelet models for turbulent combustion 

2. Effect of pressure waves on turbulent combustion 

Acoustic waves play an important role in non-reacting as well as reacting 
flows. As most experiments are performed in ducts, the acoustic modes of these 
ducts are first choice candidates to excite the flow to study. Therefore, many 
experimental shear layers are controlled mainly by acoustic excitations (Ho and 
Huerre 1984). We will start this work by looking at the effects of acoustics in 
the non-reacting mixing layer and afterwards proceed to the reacting case. In 
these two situations, the questions addressed will be essentially the same: 

- How should excitation be introduced in the direct simulation to correspond 
to experiments? Exciting the shear layer at its most amplified hydrodynamic fre- 
quency is clearly not compatible with all experimental data: many experiments 
exhibit a sensitivity to acoustics that must be considered (Masutani et al 1986, 
Poinsot et al 1987). Exciting the flow with an acoustic frequency of the duct 
might be an interesting alternative, but the best technique would be to introduce 
initial perturbations on different oscillation modes (acoustic and random modes) 
and let the system evolve without other excitation (Poinsot and Candel 1988). 
The use of reflecting boundary conditions for acoustic waves should insure that 
the flow keeps being excited by its own oscillations. It is interesting to note that 
this method has already been used successfully to predict coalescence in free 
shear layers (Grinstein et al 1987) and theoretically confirm the feedback law 
proposed by Laufer (see Ho and Huerre 1984). 

- What acoustic modes are involved in the locking phenomena? 

- How does the flow respond to acoustic excitations, and what kind of locking 
can be isolated from the results? For this study, the reacting case might be 
simpler than the non-reacting case because of generation of acoustic disturbances 
by non steady combustion in vortices. 

Many experiments have been conducted on shear layers and could be used to 
validate the results obtained from direct simulation. After a preliminary study, 
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however, the experiment of Masutani and Bowman (1986) on a reacting shear 
layer appears as a good choice because it was performed in the premixing domain 
where three dimensional effects are reduced and for Reynolds numbers which are 
accessible to direct simulation. All computations will be two dimensional. One 
of the streams will contain nitrogen and nitric oxide while the other contains 
nitrogen and ozone. Either reactant (nitric oxyde or ozone) can be injected in 
the high speed side. The Reynolds number based on the vorticity thickness and 
the velocity difference will vary between 300 and 2000. The chemistry can be 
represented by only one reaction namely: 

NO -t- O3 ^ NO^ + O 2 

This paper provides many experimental data such as mean profiles, RMS profiles, 
PDFs and power spectra density at different locations for reacting and non 
reacting cases. These data should constitute a good and critical basis for code 
testing. 

3. Validation of flamelet models for turbulent combustion 

Flamelet models for turbulent combustion are the subject of growing inter- 
est. In these models, turbulent combustion is modelled as a collection of small 
laminar flamelets which are convected, stretched, and quenched by turbulence 
without losing their laminar structure. These assumptions allow a computation 
of the mean turbulent reaction rate using only two informations: 

1- the laminar flame speed of the flamelets 

2- the topology of the flamelets ensemble 

Computing the laminar speed and the extinction limits of stretched flames has 
become a relatively easy task (Giovangigli and Smooke 1987) and different data 
bases giving flame speeds for various flow parameters are or will be available in 
the near future. Defining the topology of the flamelets is a much more difficult 
task and constitutes the basic problem of flamelet modelling. 

A powerful flamelet model has been proposed by Marble and Broadwell (1977) 
and developed by S. Candel and coworkers (1988) (see also Darabiha et al 1987). 
In this model called the Coherent Flame Model (CFM) the topology of the 
flamelets is described only by their surface E, the flame front surface per unit 
mass of gas (the flame front surface per unit volume of gas pH is more con- 
venient to use). This approach requires the chemical time to be smaller than 
the turbulent times so that the flame is essentially wrinkled and stretched by 
turbulence but that it remains a continuous interface separating burnt and fresh 
gases. An equation is written for this flame surface and coupled with the flow 
properties: turbulence field, chemical data, etc. The laminar flame speed is ob- 
tained from a data base containing computations of strained laminar flames at 
stagnation point with complete chemistry. The interesting feature in the model 
is its ability to handle chemical problems on one side to obtain the local flame 
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speed of the flamelets and the turbulence-combustion interaction on the other 
side through the flame surface equation. Although the CF model has already 
been successfully applied in different practical combustion systems, its complete 
validation will certainly be easier if direct simulation is used to ascertain most 
of the assumptions made to obtain the flame surface equation. This equation 
describes the evolution of the flame surface in each computational cell. Like the 
equations derived in turbulence modelling for the turbulent kinetic energy or the 
dissipation rate, this equation requires modelling of many physical mechanisms 
(Marble and Broadwell 1977, Candel et al 1988): 

- convection of the flamelets by the mean flow 

- diffusion of the flamelets by turbulence 

- flamelet stretching by the turbulent field (which is a source term for the 
flame surface equation) 

- annihilation of flamelets by mutual interactions (which is a sink term for the 
flame surface equation). 

+ V = EpE + Diffusion Term — Annihilation Term 
at ox ay 

The term EpE represents the flame surface increase due to the stretching of 
the flame front by the turbulent flow and is the dominant term in the flam e 
surface equation. E is the stretch rate and might be easily obtained as y/efu if 
chemistry was infinitely fast and the flame was an infinitely thin interface (e is 
the turbulent dissipation rate and v is the kinematic viscosity). This is not the 
case, and because of chemical effects, most flamelets will be quenched if they 
are submitted to high stretch rates. Therefore, one of the first difficulties is to 
give a proper estimate of E taking into account all possible interactions between 
the flame front and the turbulent field. Assumptions are usually made for each 
of these terms, but very few validations of these hypotheses exist. For example, 
Marble and Broadwell (1977) estimated E in a reacting shear layer to be of the 
order of ^ which is a proper estimate of the stretch rate of the large eddies in 
the layer but does not include any chemical effects. The interaction of a laminar 
flame and a vortex has already been studied analytically (Cetegen and Sirignano 
1987) and numerically (Candel and Laverdant 1987) but recent direct simulation 
results performed on the structure of premixed flames interacting with a vortex 
show more complete and promising results (Rutland and Ferziger 1989). 

We intend to extend these studies to test and improve some of the modelling 
assumptions used in flamelet models by studying simple interactions between 
a laminar flame front and an imposed vorticity field. This will be done on a 
diffusion flame in a shear layer or on a laminar premixed flame propagating in 
a duct and submitted to a vortex. In these two cases, the initial flame structure 
will be that of a laminar flame and the objective will be to study how this 
structure interacts with the vorticity field generated by the shear layer or by 
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the isolated vortex. The flame surface will be tracked and compared to the 
computation given by the flame surface equation of the CF model using the 
same input values for strain rates. The computations will be two dimensional. 
The main parameters will be the Damkohler number based on the chemical 
time and the roll up time of the eddy. Following problems will be more precisely 
considered: 

- In which range of Damkohler numbers does the initial flame retain its laminar 
structure? (In other words can we use the flamelet concept to describe turbulent 
combustion at this Damkohler number?) 

- What are the local values of the strain rates, and is there any evidence of 
flame extinction by strain? 

- How must the CF model equation for the flame surface be modified to predict 
correctly results given by the direct simulation? How is the effective stretch rate 
E given by the direct simulation related to the flow vorticity and to the flame 
chemical time? 

- What is the effect of eddies smaller than the flame front thickness? (High 
values of Damkohler numbers). Can the flamelet concept be extended to regimes 
where turbulence does not only wrinkle and stretch the flame but also thickens 
it? 

One may note that acoustic waves are not a first order phenomenon in this part 
of the work but that variable density effects certainly are. Density variations 
will affect the flow and hence will modify the flamelet strain rate which is one 
of the major parameters in flamelet models. It is also clear that all issues cited 
above might not be completely treated in the limited time devoted to this study. 

4. Program 

This work will be done starting with the code developed by Dr. S. K. Lele at 
NASA Ames. This code is fully compressible and it will be modified and used 
in two different Versions: 

- Version 1 will have essentially the same base than the original code, will 
deal only with non-reacting flows, and will specifically address the effects of 
acoustic waves on vortex shedding and growth. It will also be used for the 
implementation of new boundary conditions and excitation methods adapted to 
acoustic treatments and different validation tests. 

- Version 2 will have chemical reactions. This will require different modifica- 
tions of the code and also some basic validation tests such as computation of 
laminar flames. 

These two versions will be developed and used together. Since September 
1st., Version 1 has been modified to provide easy pre- and postprocessing and 
is presently used to test the code performances on acoustic waves computations 
and implement more complex boundary conditions. Version 2 is being modified 
at the same time but will not be used before first validation tests of Version 1 
will be finished. 
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Turbulent transport in the solar nebula 

By K. W. THOMPSON 

1. Problem Description 

It is likely that turbulence played a major role in the evolution of the so- 
lar nebula, which is the flattened disk of dust and gas out of which our solar 
system formed. Relevant turbulent processes include the transport of angular 
momentum, mass, and heat, which were critically important to the formation 
of the solar system. This research will break new ground in the modelling of 
compressible turbulence and its effects on the evolution of the solar nebula. The 
computational techniques which have been developed should be of interest to re- 
searchers studying other astrophysical disk systems (e.g. active galactic nuclei), 
as well as turbulence modelers outside the astrophysics community. 

2. Objectives/Milestones 

The work currently being performed system is closely connected with that of 
Cabot, Hubickyj, and Pollack (hereafter CHP), who have been using the turbu- 
lent channel flow code of Kim, Moin, & Moser (1987) to investigate incompress- 
ible turbulence in the solar nebula. The turbulence simulation code developed 
for this project will ultimately replace the central elements of the channel flow 
code, while incorporating many of the useful analysis tools developed for the 
latter. 

The objectives of this project fall into two categories: 

1. To advance the understanding of the role of compressible turbulence in the 
evolution of the solar nebula by carrying out a series of numerical experiments 
to evaluate the Reynolds stress tensor, turbulent heat transfer rate, turbulent 
dissipation rate, and turbulent energy spectrum for conditions relevant to the 
solar nebula; 

2. To advance the state of the art in turbulence simulation techniques by 
exploring new computational methods which are expected to be both efficient 
and flexible. 

The field of solar nebula physics is one of much current interest, as seen, for 
example, in the contributions to Protostars and Planets II (Black & Matthews 
1985). The field of turbulence is also one of much current interest. 

2.1 Scientific Objectives 

The first objective is scientific and largely astrophysical in nature. Turbulent 
processes are believed to have been important in the evolution of the solar nebula, 
which was a rarefied disk of gas and dust, out of which the planets, asteroids, and 
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comets formed. This disk circled the Sun during and shortly after its formation. 
The nebula is in turn thought to have formed out of the contraction of a much 
larger and even more diffuse molecular cloud. The nebula’s central star, our 
Sun, also formed out of the molecular cloud matter, and it is the Sun’s gravity 
which held the nebula together and kept it from flying apart. The combination 
of the Sun’s gravitational field, the initial angular momentum of the molecular 
cloud (retained by the nebular material), and radiative cooling is believed to have 
confined the solar nebula to a thin disk, rather than a cloud. (The formation and 
evolution of protostellar disks continue to occur in active star formation regions 
in our galaxy, and this research is relevant to nebular problems in general, not 
only to our solar system.) 

Anticipated sources of turbulence include thermal convection (driven by tem- 
perature gradients), vertical shear in the angular velocity (due to the nonuniform 
deviations from the central plane orbital motion of the gas with altitude), and 
the large velocity differential between the rotating disk matter and the infalling 
(non-Keplerian) molecular cloud material. Turbulent convection is expected to 
be strongly dependent on the compressibility of the flow; hence, the study of 
compressible turbulent flow is crucial to this project. 

The effects of turbulence on the solar nebula are many. Turbulent dissipation 
causes a loss of orbital energy and, therefore, an inward flow of mass toward the 
central star. The inflow of mass must be accompanied by an outflow of specific 
angular momentum in order to conserve total angular momentum. The inward 
transport of mass feeds most of the nebular matter to the Sun, while the outward 
transport of angular momentum causes the nebula to expand. Turbulence will 
also transport heat, both radially (parallel to the central disk plane) and verti- 
cally (in the perpendicular direction). Turbulent dissipation, mentioned above, 
is expected to be a major (perhaps dominant) source of heat for the nebula. 
Turbulent dissipation converts kinetic energy to heat, which is then radiated 
away to space, reducing the total energy of the system and making the nebula 
more tightly bound to the Sun. 

The greatest obstacle to the accurate modeling of the solar nebula is the fact 
that the length scales on which viscous dissipation takes place (and on which 
the turbulent kinetic energy is turned into heat) is many orders of magnitude 
smaller than the size of the disk. Consequently, it is not currently possible to 
create a single computational model which accurately simulates both the large 
scale structure and the small scale turbulent dissipative processes. The objective 
of this project is to simulate numerically the turbulent processes on a small scale 
and obtain a parameterization of these processes which may be used in other 
attempts to model the large scale evolution of the nebula. Much of this work 
is being done by other researchers (CHP), but their efforts must be considered 
in the developments described below, as the results of their research will be 
incorporated into this model. 

Several issues need to be addressed in order to produce a realistic model for 
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the solar nebula. The most obvious of these is the Sun’s gravity, which forces 
the gas to follow nearly Keplerian velocity profiles in its rotation about the Sun. 
Radial shear is, therefore, produced in the flow, as the orbital angular velocity 
decreases with distance from the Sun. 1 

The Keplerian flow velocities are highly supersonic in the rest frame of the 
central star. A direct simulation of the flow in the star’s rest frame is impractical, 
as the flow velocities dictate unworkably small time steps. A better approach 
is to work in a coordinate system which is comoving with the average flow in 
the model volume, so that the velocities are subsonic in the model’s coordinate 
system, and the time steps are more reasonable. 

A coordinate transformation based on the work of Rogallo (1981) will be used 
to represent the radial shear in a form which permits the use of periodic boundary 
conditions in the radial direction. The Rogallo transformation eliminates the 
need to devise boundary conditions which properly advect turbulent flow in and 
out of the radial boundaries. 

A second type of shear exists due to the matter falling in from the collapsing 
molecular cloud onto the disk. The radial and angular velocity components 
of the infalling material will be quite different from those of the rotating disk. 
While the vertical velocity component of the infalling matter will be reduced by 
the accretion shock at the cloud-disk interface, it will remain nonzero, providing 
a downward directed ram pressure on the disk material. In addition, the radial 
and singular velocity components of the infalling matter will not be affected by 
the shock and will create a large velocity shear relative to the disk material. 
This is an important question, and new techniques may have to be developed to 
study it. The effects of these factors on nebular turbulence could be important. 

A subtler gravitational effect is the vertical variation of gravity within the 
disk. The disk is very thin compared to its radial dimensions; therefore, a 
fluid element a distance z above the central disk plane experiences a downward- 
directed gravitational force (due to the central star) which is proportional to 
z. This linear variation of gravity is expected to have a significant effect on 
the convective flow. Convection in the Earth’s atmosphere takes place in an 
altitude range over which the gravitational acceleration is essentially constant. 
The nebular problem has a variable acceleration with altitude, which may give 
rise to flow patterns not seen in constant gravity environments and is expected 
to affect turbulent transport. 

Another major issue is the compressibility of the fluid. The computational 
technique of Kim et al, as currently implemented by CHP assumes an incom- 
pressible fluid, but the solar nebula is strongly stratified and is, therefore, ex- 
pected to behave like a compressible gas. Thus it is necessary to use the com- 
pressible fluid equations. Most of the work on turbulent flow in the past has 

1 Note: the radial shear by itself is stable and will not directly cause turbulence. Turbulence 
must be induced by other effects such as viscous heating and vertical shear. 



110 


K. W. Thompson 


concentrated on incompressible flows. 

The channel flow code assumes flow between two parallel plates, for which the 
fluid velocities go to zero at the plates. The no-slip boundary condition makes 
sense for flow in a channel but is not valid for flow in the solar nebula. It is 
desirable to replace these no-slip wall boundary conditions with conditions which 
are representative of the solar nebula. The nonreflecting boundary conditions 
of Thompson (1987a) are being used for this purpose. These conditions allow 
material to flow across the boundaries and allow outward propagating waves to 
leave the computational domain without generating reflections. 

2.2 Improvements to Turbulence Simulation Techniques 

The second category of objectives, which may have applications outside of 
astrophysics (e.g. the aeronautical community), focuses on the improvement of 
turbulence simulation techniques in general. 

The computer code currently in use is very flexible and is able to adjust 
to a wide variety of boundary conditions and computational techniques. It is 
designed to be easily modified so that new features may be added and differ- 
ent computational techniques incorporated as needed. For example, the spa- 
tial derivative methods can be modified at will (e.g. from finite difference to 
pseudospectral) by altering only that section of the code which computes these 
derivatives. 

The techniques in use are described in section 3. 

2.3 Milestones Passed 

The work performed so far has been concerned with analyzing alternative 
computational methods, developing a suitable computing strategy, and validat- 
ing the computer simulation. Early design and testing focused on simple linear 
and nonlinear wave simulations, for which error analysis and convergence testing 
could be performed easily. The later stages have involved the full set of compress- 
ible fluid equations. Successful tests have validated the design of the computer 
code and are described below. One of these tests has had the added benefit of 
generating a significant new research project in its own right, as described below. 

The major program tests are: 

1. Rarefaction wave. 

2. Oscillating atmospheric column. 

3. Hydrodynamic escape. 

4. Unstable shear flow (Kelvin-Helmholtz instability). 

2.3.1 Rarefaction Wave 

The rarefaction wave is the well known solution to the adiabatic expansion 
of a uniform gas into a vacuum in the absence of gravity. The exact solution 
to this problem has been given in many places, (e.g. Landau & Lifshitz 1959, 
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Thompson 1986). Time dependent numerical simulations of the continuous part 
of the rarefaction wave were highly accurate and exhibited the proper fourth 
order convergence with grid refinement. 

2.3.2 Oscillating Atmospheric Column 

This problem describes a column of atmosphere in a constant gravitational 
field. The density and pressure drop off exponentially with altitude z , as e~ z ! h , 
where h is the vertical scale height. This is a one dimensional problem, as the 
solution does not depend on x or y. The initial state is isothermal and is stable. 

Just as in an organ pipe, one can set up standing waves at certain discrete fre- 
quencies. I have set up an adiabatic perturbation to excite the lowest frequency 
mode of the column for problems ranging from 1 to 50 vertical scale heights 
in extent. In all cases the time dependent solution accurately reproduced the 
expected oscillatory mode and displayed the correct fourth order convergence. 
The exponential dependence of the unperturbed state makes this problem an 
obvious choice for solution by the methods of section 3.1. 

2.3.3 Hydrodynamic Escape 

This problem is also one dimensional but spherically symmetric. It describes 
the escape of a planetary atmosphere into space. Gravity obeys the inverse 
square law, and the flow velocity goes smoothly from highly subsonic near the 
planet’s surface to highly supersonic at large distances. Steady state solutions 
have been known analytically for some time (Bondi 1952). It is an interesting 
fact that the same equations describe the solar wind and the accretion of gas 
onto a protostar embedded in a spherically symmetric cloud. 

Although the analytic solution is known, this problem presents special prob- 
lems for numerical techniques (Zahnle 1988) . I have found it to be a challenging 
problem (and a useful one for uncovering some subtle boundary condition prob- 
lems) but can now produce stable and accurate solutions. The problem contains 
enormous pressure and density gradients near the lower boundary, and the den- 
sity and pressure vary by several orders of magnitude in the model volume. This 
problem is rendered tractable by the logarithmic techniques described in section 
3.1. Previous attempts to solve this problem numerically with more standard 
methods have failed. 

The more general problem of hydrodynamic escape of planetary atmospheres 
which are affected by thermal conductivity and solar heating is of considerable 
interest to researchers in the field. The successful solution of the more basic 
problem will serve as a starting point for further research on this subject by 
Zahnle, my co-investigators, and myself. 

2.3.4 Unstable Shear Flow 

This is a well known problem in the atmospheric sciences, and one which is 
two dimensional in nature. We have an initially isothermal atmosphere whose 
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density and pressure vary exponentially with altitude z as e~ x ! h (rectangular 
coordinates). The atmosphere contains a shear flow with velocity components 
(U(z), 0,0), where U(z) = U 0 tanh(z/d), d < h, and Uq is subsonic. We perturb 
the flow field by adding small velocity perturbations u(x, z) and w(x,z) so that 
the perturbed velocity is ( U + u,0, u;). Perturbation theory (Chandrasekhar 
1961) shows that the stability of the flow is determined by the Richardson number 
Ri = gcP/hU^. A Kelvin-Helmholtz instability arises if Ri < 1/4, while the flow 
is stable if Ri >1/4. 

Calculations performed with Richardson numbers of 0.05 and 0.15 do indeed 
give rise to unstable flow. These calculations show that the initially small pertur- 
bations grow and create a vortex. Similar calculations with Richardson numbers 
of 0.35 and 1.25 are stable and do not lead to vortex formation. These results 
are consistent with the predictions of linear stability analysis. 

3. Technical Approach 

The equations to be solved are the compressible Navier Stokes equations. 
They describe the time evolution of a compressible fluid in three dimensions and 
incorporate the effects of a variable gravitational field, radial shear, molecular 
viscosity, thermal conductivity, and internal heat sources. The solution process 
is made challenging by the large variation in density and pressure which occurs in 
this problem. A new technique has been developed to handle the range problem 
as described below. 

The purpose of this calculation is to study relatively small scale turbulence 
in order to characterize the effects of turbulence on the overall flow without 
resorting to ad hoc turbulent viscosity approximations. Therefore, the actual 
volume to be simulated consists of a small “box” of nebular material in the disk 
of relatively small radial and angular extent as seen from the Sun. 2 The box 
includes one to several scale heights of pressure variation in the direction per- 
pendicular to the disk midplane. Periodic boundary conditions are used at the 
boundaries facing “sideways” into the disk material (i.e. the (f> boundaries) . The 
r boundaries make use of Rogallo’s method (Rogallo 1981) to specify periodic 
boundary conditions. The d (vertical) boundaries require non-trivial boundary 
conditions and are currently handled by the nonreflecting conditions described 
by Thompson (1987a). 

3.1 The Range Problem 

The density and pressure in the nebula fall off very rapidly with altitude above 
the nebula midplane, roughly as , where h is a scale height. 3 Thus the 

density and pressure may vary by orders of magnitude throughout the model 
volume, which poses a difficult challenge to conventional numerical techniques. 

The spherical coordinate system is assumed for this discussion. 

o 

The unusually rapid falloff is due to the variable gravity, which increases with altitude. 
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I have developed a technique to deal with scalar fields which vary over such a 
large range by using the logarithms of density and pressure instead of the basic 
quantities themselves. For example, the one dimensional continuity equation, 


dp , d 


is replaced by 


dlnp 3 In p du 

HT + u ~aT + ai = 0. 


and similarly for the energy equation. Thus an exponential dependence of p 
on x is reduced to a polynomial dependence of In p on x, which the numerical 
scheme can cope with more readily. All derivatives of density, pressure, and 
thermal energy may be replaced by their logarithmic forms as above, leading 
to a more tractable system. An added benefit is that the density and pressure 
cannot become negative with this approach. 


3.1.1 The Numerical Techniques 

The current code uses standard fourth order finite difference formulas to eval- 
uate spatial derivatives in all directions (see, for example, Thompson 1987b). 4 

Time integration is handled by the usual fourth order explicit Runge-Kutta 
method, which is simple to implement and has excellent stability properties. 
This combination of spatial derivative and time integration methods provides 
global fourth order convergence. The fourth order convergence rate has been 
verified repeatedly on a large number of test runs and implies that this code 
requires significantly fewer grid points than a second order code to achieve a 
given level of accuracy. 

The choice of an explicit method over an implicit method stems from the need 
to resolve the smallest features present in the flow. At the smallest length scales, 
viscosity dominates the evolution of the flow. Since we need to simulate the 
dissipation of kinetic energy to heat accurately at these scales, the grid spacing 
and time steps necessary are set by the properties of the flow and are the same 
whether explicit or implicit methods are used. The optimal grid spacing and 
time step are those for which the propagation and viscous Courant numbers 
are equal. Consequently, the simpler (and faster) explicit approach has been 
selected. 


4 Compact difference formulas of the compact (Pade) type (Rubin & Koshla 1977) have been 
evaluated but have not been found to offer significant improvements over the standard fourth 
order formulas. 
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Transition to turbulence 
in laminar hypersonic flow 

By J. VAN DER VEGT 


1. Introduction 

This report gives a short discussion of the progress in a recently started project 
aimed at the prediction of transition to turbulence in hypersonic flow. The pre- 
diction of transition to turbulence is a very important issue in the design of space 
vessels. Two space vehicles currently under investigation, namely the aeroas- 
sisted transfer vehicle (AOTV) and the trans-atmospheric vehicle (TAV), suffer 
from strong aerodynamic heating. This heating is strongly influenced by the 
boundary layer structure. These aerospace vehicles fly in the upper atmospheric 
layer at a Mach number between 10 and 30 at very low atmospheric pressures. 
At very high altitudes the flow is laminar, but when the space vessel returns 
to a lower orbit, the flow becomes turbulent and the heating is dramatically 
increased. 

The prediction of this transition process is commonly done by means of ex- 
periments. The experimental facilities available nowadays cannot model the 
hypersonic flow field accurately enough by limitations in Mach and Reynolds 
number. These facilities also have a large free stream disturbance level which 
makes it very difficult to investigate transition accurately. An alternative ap- 
proach is to study transition by theoretical means. Up to now numerical studies 
of hypersonic flow only discussed steady laminar or turbulent flow. The project 
discussed in this report tries to extend this theoretical approach to the study of 
transition in hypersonic flow by means of direct numerical simulations and addi- 
tional theoretical investigations to explain the mechanisms leading to transition. 
This report gives a brief outline of how we plan to carry out this research. 

2. Important Physical Phenomena 

There are several important new physical phenomena in hypersonic flow which 
are not important in supersonic flow. The flow exhibits extreme heating, espe- 
cially in the stagnation point region. This causes non-equilibrium chemical re- 
actions, dissociation, and ionization, which have to be taken into account. Also, 
the effect of radiation cannot be neglected. The gas is rather dilute, so the lim- 
its of the continuum approach are reached, but with slight modifications, the 
Navier-Stokes equations still can be used. 

3. Transition to Turbulence 

There are several techniques for studying the transition phenomena. The 
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most widely used technique is linear stability theory, using the normal mode ap- 
proach in both space and time coordinates. This technique is, however, not very 
reliable for the prediction of the position and behavior of the actual transition 
phenomena and, in many cases, additional empirical information is necessary 
for accurate prediction. This empirical information is very scarce in hypersonic 
flow so accurate numerical solutions are necessary. 

The main reason for the failure of linear stability theory can be attributed to 
the fact that finite amplitude disturbances, requiring non-linear stability theory, 
cannot be neglected. This is demonstrated by experiments on a cone conducted 
by Stetson (1988). A more general valid theory describing the transition process 
is given by Ruelle and Takens (1971). It assumes that the change from a lami- 
nar to a turbulent state takes place by means of a finite number of bifurcations. 
Studying these bifurcations by theoretical means and direct numerical simula- 
tions can give more detailed information about the transition process. This 
method allows large disturbances and gives more freedom than a linear theory. 

Direct numerical simulations of hypersonic flow accurate enough to predict 
transition will be very time consuming. An additional problem with direct nu- 
merical simulations is caused by the fact that there axe so many independent 
parameters that it is hardly possible to investigate all relevant situations by this 
method. It is, therefore, necessary to study the transition problem by means 
of a coupled theoretical and numerical approach. The aim of this study is to 
obtain detailed information on the actual mechanisms leading to transition in 
hypersonic flows and to obtain information useful for engineers. 

4. Achievements 

The project started in September 1988, so there are only limited results. Dur- 
ing this period, the literature and basic phenomena in hypersonic flow were 
studied. The results of this study are to be used to make a project proposal for 
the study of hypersonic transition. 

5. Proposal for the Calculation of Transition to Turbulence in Lami- 
nar Hypersonic Flow 

General: The main purpose of this research is to get insight in the phenomena 
which lead to transition from laminar to turbulent hypersonic wall bounded 
flows. In order to accomplish this goal, two main tools will have to be developed: 

- An accurate code for the direct numerical simulation of time dependent hyper- 
sonic flows. 

- A mathematical model for investigating the stability and bifurcation phenom- 
ena leading to transition in hypersonic flow. 
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Due to the facts that we have to start from scratch and that there are large 
uncertainties about the method, it is very important to follow a stepwise proce- 
dure. This approach is summarized below: 

-A Solution of Two-Dimensional Navier-Stokes Equations for a Flat 
Plate 

The first step will be to make a computer model for the viscous supersonic 
two-dimensional flow on a flat plate at low Mach number (Af < 3) using a finite 
difference method. 

It is assumed that there are no shocks in the flow field. The numerical model, 
however, must solve the Navier-Stokes equations in conservation form together 
with the ideal gas equation of state. The use of a conservative formulation 
provides the opportunity of shock capturing in a later stage of the project. The 
numerical method has to be at least partially implicit in time in order to be 
able to add a chemistry model for the reacting gas flow while avoiding serious 
restrictions on the time step due to stiffness. Possible numerical methods are 
the implicit Me Cormack and Beam and Warming schemes which are second 
order accurate in time and space (Anderson, Tannehill and Pletcher 1984). 

Due to the geometric simplicity a new program will be written instead of using 
an existing general purpose computer program which would need modifications 
to be efficient for a flat plate. 

• The purpose of this part is: 

1. Obtain experience with these type of calculations. Although existing numer- 
ical methods will be used, it is important to test these methods on a relatively 
simple problem and especially to investigate the effect of numerical parameters 
like time step, number of grid points, stability, and numerical diffusion. 

2. Test the accuracy of the numerical algorithm by comparing the growth rate 
of an Orr-Sommerfeld wave with the prediction obtained with linear stability 
theory by Mack (1984). 

3. The direct numerical simulations will be used to compute the evolution of 
finite amplitude disturbances and study the effect of wall temperature. 

-B Three-Dimensional Direct Numerical Simulations of the Flow on 
a Flat Plate 

There is a severe restriction on using a two-dimensional model for direct nu- 
merical simulations because, at high Mach numbers, oblique waves become more 



118 


J. van der Vegt 


unstable than two-dimensional waves. In order to investigate this phenomenon, 
the numerical model will be extended to three-dimensional flow on a flat plate. 

• This part of the research should give a more accurate representation of the 
flow field and offer the opportunity to investigate three-dimensional effects. It 
can also be used for a thorough investigation of the accuracy of the e N method 
widely used in engineering applications. 

-C The Effects of Chemistry 

The effects of chemistry on transition are not very clear. In order to investigate 
these phenomena, both with respect to their effects on the numerical algorithm 
and on the flow field, they will be added to the two-dimensional numerical model. 
Initially the chemical model for non-equilibrium chemistry of Candler (1988) will 
be used. This model is a seven species chemical model for N 2 , O 2 , NO, N, O, 
NO + and e~. The effects of ionization will be neglected in this stage, so the 
rate equations for NO + and e~ will not be used. 

• This part of the research should give information about the effects of the 
chemical reactions on transition phenomena and velocity profiles. 

-D Further Study of Transition 

In the investigation of the phenomena leading to transition, linear stability the- 
ory is not always accurate enough. The most important reason is the fact that 
linear stability theory gives no information about finite amplitude disturbances. 
A non-linear stability study will be conducted to investigate these effects. 

• This part of the study would give more detailed information about the actual 
phenomena leading to transition and they can be further investigated with the 
numerical model. 

-E The Effects of Internal Energy and Ionization 

At higher temperatures, energy associated with internal degrees of freedom 
and ionization effects become important. These effects are incorporated by 
adding an additional energy equation for the vibrational energy and an equation 
for the electron translational energy. At this stage, it is supposed that there 
are three relevant temperatures, viz. a translational temperature, a vibrational 
temperature and electron temperature. 

• This part of the research should give information about the internal energy 
and ionization effects on transition. 
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-F Direct Numerical Simulation of the Flow Field around a Cone 

The procedure discussed for the flat plate can also be used to investigate the 
flow field around a cone. The calculation of the flow field around a cone is 
more difficult, so it is advantageous to start this research when more experience 
is obtained with the flat plate. An interesting phenomenon in this flow is the 
interaction between the boundary layer and a shock. The shock has a strong 
effect on the flow field. However, it requires a high resolution which is probably 
difficult to obtain in three-dimensional flows. The study will, therefore, start 
with an axisymmetric model. 

• This part of the research can provide more information about the effects of 
interaction between shock and boundary layer on transition. 
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Direct numerical simulations of 
turbulent convection with a variable 
gravity and Keplerian rotation 

By W. CABOT 


1. Motivation 

Thermal convection has been proposed as a possible mechanism for genera- 
tion and maintenance of turbulence in the inner accretion disk regime of the 
primordial solar nebula. Resulting Reynolds stresses would produce a torque on 
the disk, causing “viscous” spreading and the eventual dispersal of the gas com- 
ponent of the disk; this process was coeval with and coupled to the formation 
of the planetesimals and, ultimately, planets. Turbulent convection under solar 
nebula conditions has been described to date by ad hoc and mostly untested 
models of thermal convection and turbulence. The models of Lin & Papaloizou 
(1980) and Cabot et al (1987), in particular, give vastly different results in terms 
of the vertical distribution and intensity of turbulence, the efficiency of convec- 
tion, and structural stability properties. There are yet few, if any, conclusive 
astronomical observations of solar nebula analogues, nor does indirect evidence 
from planetary compositions provide clear constraints on these models. It is, 
therefore, of fundamental interest to design experiments with the basic physical 
features of the solar nebula in order to constrain old and new models. Although 
solar nebula conditions cannot be reproduced in the laboratory, numerical sim- 
ulations of hydrodynamic flows, which have been very successful in describing 
aerodynamic flows, can be suitably modified to provide “experimental” data for 
solar nebula modelling. 

2. Objectives 

The goals of this project are (l) to modify an extant, “proven” hydrodynamics 
code with the most important features of the solar nebula and other thin accre- 
tion disks: buoyancy terms to generate convection, internal heating representing 
the release of gravitational potential energy, a variable gravity linearly propor- 
tional to the distance from the vertical midplane due to centrifugal balance, 
rapid rotation with axis aligned with gravity, and Keplerian rotational shear; 
(2) to determine the effect that these features have on the turbulent convec- 
tion by introducing them individually and to determine the cumulative nature 
of the turbulent convection for accretion disk conditions; and (3) to model the 
convection (viz., the convective heat flux) and the turbulence (viz., turbulence 
intensities, heat dissipation, and Reynolds stresses). In this manner, prior solar 
nebula models can be tested and their deficiencies rectified. 
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3. Work to Date 

Uniform rotation. The code for direct numerical simulations of a incompress- 
ible flow in a semi-infinite channel (Kim, Moin & Moser, 1987) was modified 
to include buoyancy in the Navier-Stokes equation in the Boussinesq approxi- 
mation, the variable gravity, and uniform rotation. Both no-slip and no-stress 
boundary conditions on the vertical walls were implemented (the latter being 
more realistic for the solar nebula). Sequences of simulations with varying ro- 
tation rates were performed, from no rotation to rates approaching the limit 
of stability. A linear analysis was performed to determine theoretical critical 
rotation rates for marginal stability. It was found that the convective fluxes 
in the numerical simulations become negligible between the marginal limits of 
stationary and oscillatory convection, and there may also be finite amplitude 
effects occurring in this rotation regime. 

The turbulence intensities and efficiency of the convection were found to be 
virtually undiminished in the midchannel region where buoyancy vanishes; this 
region posed a major uncertainty in prior solar nebula modelling, and this result 
not only clarifies the nature of the turbulence there but allows some simplifi- 
cation in modelling. An extension of the standard stellar mixing length model 
for convective heat fluxes and speeds was developed to include viscosity and 
rotation and was found to agree qualitatively with the rotational dependence of 
the simulation results, although convective fluxes and speeds in the model were 
found to be a factor of 3 too low for standard values of coefficients. Canuto & 
Goldman’s (1985) turbulence model, used in Cabot et aV s (1987) solar nebula 
model, also gave a qualitatively similar dependence of convective heat flux and 
speed on rotation, but gave convective fluxes that were SO times too low. 

The vertical distribution of convective buoyancy production of turbulence ki- 
netic energy has a bimodal shape with peaks on either side of midchannel and 
vanishing at midchannel. For low rotation rates, nonlinear diffusion (primarily 
by convective transport) redistributes the turbulence so that it is dissipated with 
a flat interior profile. For high rotation rates, the overall diffusion is suppressed 
(with diffusion by pressure fluctuations now dominating) such that dissipation 
near midchannel is depressed. This result further supports the assumption made 
by Cabot et al (1987) that heat dissipation is evenly distributed throughout the 
convective region, than it does the model by Lin <fc Papaloizou (1980), in which 
the turbulence intensities and dissipation is sharply peaked at the midplane. 
Since the heat dissipation is ultimately the heat source in the solar nebula, it 
is important to assess the sensitivity of the heat dissipation distribution to that 
of the internal heat source. Tests with centrally concentrated internal heating 
show that the interior heat dissipation profile is very insensitive to that of the 
heat source, which is again consistent with the assumption that the dissipation 
will have a relatively flat interior distribution. 



Turbulent convection with variable gravity and Keplerian rotation 


123 


Rotational shear. Linear rotational shear is treated in the numerical integra- 
tion by working in a cosheared reference frame, in which the governing equations 
remain horizontally homogeneous, and implementing Rogallo’s (1981) remesh- 
ing transformation in the homogeneous horizontal directions, which allows one 
to circumvent the tendency of the numerical grid to become overly distorted as 
it follows the sheared flow. The remeshing transformation was implemented in 
the channel code and has been extensively tested. The tendency of the stream- 
wise length scales to become elongated and the spanwise length scales (in the 
direction of the shear) to become compressed requires balancing horizontal box 
and mesh sizes in the simulations in order to optimize the numerical resolution, 
and runs are currently under way in order to determine the optimal numerical 
domain(s). Several simulations were carried out with less than optimal but ac- 
ceptable resolution, which have provided the preliminary results subsequently 
quoted. 

Linear stability analysis predicts that the rotation profile is centrifugally sta- 
ble if the specific angular momentum gradient increases with distance from the 
rotation axis (5 > — 2H, where S is the shear rate and Cl the rotation rate). 
This limit was tested and verified with the numerical code with no buoyancy 
forces in effect. This indicates that the rotational shear flow for Keplerian ro- 
tation (5 = — 3fl/2) is centrifugally stable and, by itself, cannot maintain the 
turbulence in accretion disks. Another verification of this proposition was found 
by allowing a convecting, rotationally sheared flow simulation with well devel- 
oped Reynolds stresses due to the shear to decay when the internal heating was 
extinguished; both the convective and shear production rates decayed to zero, 
indicating that it is convection alone that drives the turbulence in this flow and 
indirectly maintains Reynolds stresses by the “catalytic” effect of the rotational 
shear. 

Preliminary results for Keplerian rotation at a moderate rotation rate with 
both no-slip and no-stress (coshearing) vertical walls show well developed 
Reynolds stresses in the horizontal components (— uuJ). The corresponding shear 
correlation coefficient was found to be 0.21 for no-stress walls and 0.23 for no- 
slip walls, compared to 0.45 for homogenous, unidirectional, plane shear flows 
(Townsend, 1956); these values were very uniform across the channel, even at 
the walls. The shear production rate of kinetic energy was about 16% of the 
convective production rate. Simulations at different rotation rates are needed to 
determine the robustness of these results. For uniform rotation with the same 
inertial frequency, the rms vertical velocity (v rms ) was found to exceed the equal 
streamwise (u rm „) and spanwise (uVm») components by about 1.4, because it is 
the component directly receiving the production by convective buoyancy. For 
the Keplerian rotational shear case, v rmt « 1.3u; rma and w rm , « 1.4u rmg , 
compared to the result for unidirectional, plane shear flow of u rmt > v rm , > 
u) rma . Because convective buoyancy is driving the turbulence in the numerical 
simulation, one expects v rms to still be the largest. The discrepancy in the 
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horizontal components may be due to rotation: The production rate in the 
streamwise component is — (2fi + S) uw — — fi/2 tluJ < 0, rather than —Suw > 0 
without rotation, and the spanwise production rate is 20 uw > 0. 


4. Future Work 

In upcoming work, numerical simulations using sequences of rotation rates 
with Keplerian shear will be constructed. The results will be used to test the 
robustness of the value of the shear correlation coefficient and the ratio of shear 
to convective production, or to suggest ways to model whatever variations arise. 
The stability behavior of the rotationally sheared flow at very rapid rotation 
rates will be compared with uniform rotation results to determine any effects of 
the shear; also the effect of shear on the modelling of convective fluxes and speeds 
will be examined. The interior distribution of kinetic energy dissipation with 
both convective and shear production will be examined for different rotation 
rates; modelling for the dissipation and the Reynolds stress will be attempted. 
Simulations with different Rayleigh and Prandtl numbers will be performed in 
order to assess their effects on the preceding key quantities with the viscosity- 
independent product of Rayleigh number and Prandtl number held close to solar 
nebula model values (~ 10 4 to 10 s ) ; lower Prandtl number simulations will allow 
higher, more realistic rotation rates to be accessed. At this stage, much more 
physical and realistic models and paxametrizations of turbulent convection in 
thin accretion disk environments, such as the solar nebula, should be attainable 
than have been available heretofore. 

In the realm of the more distant future are a number of modifications that 
would make the numerical simulations more realistic. The problem of unrealistic, 
impermeable vertical boundaries and the imposed vertical scale of the convec- 
tive region can possibly be relaxed by introducing compressibility to the code, 
at least in the form of a variable vertical density. For coverage of many density 
scale heights, the vertical boundary region would be buffered by a convective sta- 
ble region and the convective motions would select their natural vertical scale; 
however, implementation of this change will require extensive modifications of 
vertical integration scheme in the present code, probably requiring finite dif- 
ferencing rather than the present Chebyshev polynomial expansion. Since the 
dissipation of kinetic energy should properly be equated with the internal heat 
source, ways to make a self-consistent “bootstrap” will be explored; this will 
provide a self-consistent way to specify the energy source, and will indicate if 
the coupling between energy source and turbulence leads to any unanticipated 
behavior. 

This work is being carried out in collaboration with the Space Science division 
at NAS A/ Ames (J. Pollack and P.Cassen) and with NASA/GISS in New York 
(O.Hubickyj and Y.Canuto). 
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Experimental investigation of flow over 
a backward facing step-progress report 

By L. W. B. BROWNE 1 


Introduction 

This report is a brief statement of the work contributed by the author to the 
backward-facing step experimental facility located in Test Cell No. 3 in the FML 
Building, N-260, at the NASA-Ames Research Center, during his sojourn with 
the Center of Turbulence Research — 3rd February to 28th August 1988. The 
work was carried out in close cooperation with Dr. Srba Jovid, but since this is 
not a formal report of the Center of NASA-Ames, it has only a single author. 

1. Wind Tunnel 

Figure 1 shows the sizes and arrangement of the wind tunnel used for these 
experiments. 

The sizes shown are for the final arrangement of the tunnel. They give an 
area contraction ratio of 9.4, a tunnel aspect ratio (width/height) of 2.17 and a 
step aspect ratio (tunnel width/step height) of 11.0. 

In order to use temperature as a passive marker of the flow, it was considered 
desirable to have the floor of the tunnel heated to a uniform temperature of about 
10° C above the free-stream temperature. This was accomplished by placing 
6.3mm thick aluminium plates, with heaters, on the tunnel floor. These rest 
on spacers placed on the original plexi-glass floor and occupy the 1690mm and 
2310mm lengths shown in Figure 1. The aluminium sheets have thin-foil heaters 
glued to their underside — the heaters being arranged to cover the whole surface 
of the aluminium. Four heating zones were provided, the voltage to each zone 
being controlled manually using a variac. 

Approximately 40v A.C. has been found to give the required 10° C temperature 
difference when the free stream velocity is lO.Om/s. The plate temperatures are 
monitored using 13 thermocouples located on the center-line of the plate in small 
holes drilled to a depth of 5mm. A temperature distribution with a variation of 
less than *0.2°C over the whole plate can be achieved. 

To provide a smooth entrance from the contraction section to the aluminium 
surface, the contraction section was raised until the aluminium surface and the 
contraction surface coincided. The resulting gap at the top of the tunnel was 
filled with a free-formed aluminium sheet held in place with duct tape. 

1 Permanent Address: University of Newcastle 
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FIGURE 1. Layout and Sizes of Wind Tunnel 

The flow was tripped, at the start of the aluminium plate, using a 1.6mm 
diameter rod followed by a 110mm width of 40 grit emery paper. 

To provide an approximate zero pressure gradient flow, the sides of the tunnel 
were angled slightly outwards. 

2. Air Flow 

After the difficulties associated with trying to use the large (240,000 c.f.m.) 
FML compressor to provide a suction flow through a sonic nozzle, it was decided 
to install a standard suction fan system. 

A general purpose fan, manufactured by the New York Blower Company, was 
purchased. This has a 460mm diameter fan wheel with airfoil-shaped blades 
capable of running at 2100 rpm. At a speed of about 800 rpm the required tunnel 
velocity of 10.0 m/s was achieved. As discussed later, this velocity provides an 
R.0 for the flow approaching the step of 1680. The fan is driven by an “Adjusto- 
Spede” eddy current clutch unit of the Eaton Corp. The motor (2 hp) runs at 
a constant 1130 rpm while a controlled current supply to the clutch provides 
variable speed at the output shaft. A model 3000 controller from the same 
company controlled the current. 

Only 1/2 hp is required for the fan under our operating conditions, but the 
larger unit will allow other operating conditions to be considered. The single- 
turn potentiometer of the controller was replaced with a 10-turn unit with 
counter and was located adjacent to the experiment control area. This pro- 
vides excellent speed control — the free stream velocity can be held steady, 
within + 0.2%, over typically several hours of operation. 

A transition piece between the end of the tunnel and the inlet to the fan was 
built (see Figure l). 

3. Cold- Wire Anemometers 

Constant-current anemometers are commonly used with fine wires (typically 
0.63/im diameter) to measure temperature fluctuations in turbulent flows. Ian 
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FIGURE 2. Circuit diagram of constant-current anemometer. (Designed by 
I. Miller — University of Newcastle, Australia) 

Miller of the Dept, of Mechanical Engineering at the University of Newcastle 
has evolved, over the years, an excellent design for such units — See Figure 2. 

Twelve circuit boards, one complete unit (to use as a guide), circuit diagrams, 
board layouts and instructions were supplied from Newcastle by Ian Miller and 
subsequently twelve constant-current anemometers were built at NASA-Ames. 
Some time was spent correcting problems (faulty components) but subsequent 
trials and use have shown the units to be operating very satisfactorily. 

4. Hot-Wire Anemometers 

Constant-temperature anemometers have also been developed by Ian Miller 
at the University of Newcastle — see Figure 3. 

The signal to noise level of these units is as good as the best commercial units 
and their high frequency response (« 25kHz) is quite adequate, although not as 
high as can be achieved with the commercial units. Three complete units, five 
circuit boards, circuit diagrams, etc. were sent from Newcastle by Ian Miller. 
The three units are working very well, but due to parts supply problem, the 
other five units have not been completed. It may take another month to finish 
the five units. When these have been tested satisfactorily, the original three 
units will be returned to Newcastle. 

5. Cold- Wire Rakes 

Rakes of cold wires often provide the best means for the detection of coherent 
structures in a turbulent flow in which temperature has been added as a passive 
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FIGURE 3. Circuit diagram of constant-temperature anemometer. (Designed 
by I. Miller — University of Newcastle, Australia) 

scalar. Units, built at the University of Newcastle, and twelve plug-in probes 
were brought to NASA-Ames by the author. The units allow a 2-D array of 
wires to be used if required. 

6. Wire Preparation 

A facility for cold and hot wire preparation (soldering and etching) was de- 
veloped. 

7. Miniature 3-Wire Probe 

To investigate near-wall structures, both in the upstream-of-step region and 
the downstream-of- reattachment region, it was felt that a miniature x-probe 
with an associated cold wire would be useful. Such a unit was built. The 
hot wires are 0.3mm long (« 3 Z* where L* is the Kolmogorov length scale) 
separated by 0.2mm. The wires, 2.3/i m diameter platinum, are etched from 
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Wollaston wire and attached directly to the prongs. Their Lfd of 130 is a little 
small, but from other work, were considered to be satisfactory. 

To allow for more accurate measurements in the high turbulence intensity 
regions of the flow, the included angle between the wires was made 140° rather 
than the usual 90°. The cold wire is 0 . 63 /im Pt - 10% Rh etched from a 1mm 
length of Wolluston wire soldered to prongs that can be moved forward and away 
from the x-wires to provide for easiest replacement of the wires. In its operating 
position, the cold wire is located in front and at right angles to the hot wires 
and as close as possible to them. The etched length of 0.5mm does not interfere 
with the flow over the hot wires. 

8. Software 

A number of programs were developed. These relate to instrument calibration, 
data taking, and data processing. 

9. Results 

9.1. Mean Free-Stream Velocities 

Using a Pitot-tube it was found that the free-stream velocity, of lOm/s, out- 
side the boundary layers on all surfaces of the tunnel, was uniform at any x 
position with *0.2%. There was a slightly favorable pressure gradient from the 
trip to the step, the center-line-mean velocity increasing by 0.3m/s over that 
distance (1690mm). A reference point for measuring the free stream velocity 
was established at 300mm upstream of the step. 

9.2. Spanwise Velocity Variations in Boundary Layer Upstream of Step 

The spanwise velocity variations in the floor boundary layer were typical of 
wind tunnel flows. Figure 4 shows the velocity variations, obtained with a Pitot 
tube, at x = —225mm (i.e. 225mm upstream from the step), at two distances 
from the wall. 

The same patterns exist right through the boundary layer at any x position 
between the trip and the step (see Figure 5). 

Note that the peaks and valleys of the variations occur consistently at the 
same z locations at all x positions and at any position in the boundary layer. It 
is also interesting to note that the same velocity variation patterns exist in the 
boundary layer before it arrives at the trip (see Figure 6). 

Swearingen and Blackwelder (1987) found similar patterns of spanwise varia- 
tions of the mean velocity in a wind tunnel boundary layer. 

Bradshaw (1965) and Mehta and Hoffmann (1987) carried out investigations of 
this phenomena. The results they obtained changed whenever the last screen in 
the screen box was changed (or even when the last screen was vacuum cleaned) , 
but the spanwise variations were always present. The best results obtained by 
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FIGURE 4. Spanwise distribution of mean velocities at x = —225mm at two 
positions in the boundary layer. 


Mehta and Hoffmann (1987) was a + 5% peak to peak variation in the spanwise 
values of Cj, although a +17% variation with earliest screens had been mea- 
sured. We obtained a +7% variation in the spanwise values of Cj at x/H « —6 
(H is the step height), see Figure 7, and considered this to be satisfactory. 

It was decided to carry out all future measurements at a z of 252mm. This 
point was still near the centre-line position ( z — 208mm) and lay on the Cj 
curve at a point where the values were not changing too rapidly. 

Note that (Mehta, verbal communication), since the spanwise variations are 
due to the final screen, (essentially there is some agglomeration of the streamwise 
vortices formed by the screen) a boundary layer developed on a plate placed 
centrally in the wind tunnel still has similar spanwise variations in the velocity 
and Cj values. 

Recently, at a Berlin Conference (Second European Turbulence Conference, 
Aug. 30-Sept. 2, 1988, Techmische Universitat Berlin) Dengel and Fernholz 
(1988) reported that they had reduced Cj variations, caused by screens, from 
+ 15% to 1% by using special “well behaved and carefully selected screens.” 
Their final screen was produced by stamping rectangular holes in a piece of sheet 
metal — the sheet being about 0.5mm thick and the holes having a size of about 
8mm x 5mm with their centre-lines being on a grid of about 11mm x 8mm. 
The solidity was thus about 0.4 - close to the norm for wind tunnel design. 
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FIGURE 5. Spanwise distribution of mean velocities 355mm downstream of 
the trip at three positions in the boundary layer. 



FIGURE 6. Spanwise distribution of mean velocities 10mm upstream of the 
trip at one position in the boundary layer. 

9.3. Fully-Developed Nature of the Flow Upstream of the Step 

In order to determine whether the flow had become “fully developed” as it 
approached the step, and also to check some of the instrumentation, it was 
decided to obtain profiles through the boundary layer at x/H « —6 (actual x 
position was -225mm) using a Pitot tube, a single hot wire, the miniature X- 
probe and a standard Disa X-probe (wire d — 5pm, l — 1.25mm). To ensure 
that there was no upstream influence of the step at this point, the step was later 
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FIGURE 8 . Mean velocity distribution across the boundary layer at i /H & —6. 
• Pitot; x Single hot wire; o Miniature X-probe; A Disa X-probe. 

“removed” by placing a length of ply-wood, in the downstream section of the 
tunnel, level with the upstream section, and the profile measurements were then 
repeated. 

All mean velocity U profiles, plotted using outer layer values for non- 
dimensionalizing, gave good agreement — Figure 8. 

To compare the velocity profiles using wall variables, it was necessary to de- 
termine the local Cj value. This was obtained by three different methods: 
a) Using a Preston tube. The one we used was a flattened tube of height 
0.252mm. (Subsequently, it was found that this corresponded to a height of 
y + « 7). The data was reduced using the Patel (1965) correlation equations. 





Flow over a backward facing step 


135 



FIGURE 9. Mean velocity distribution across the boundary layer at x/H « —6. 
• Miniature X-probe with step; o Miniature X-probe without step; A Single 
hot wire; x DIAS A X-probe; — Best fit in log region. 

b) Using the Clauser (1954) charts. These are based on 

v * = + 51 

c) Using the Coles (1962) curves of R$ vs Cj. From the mean velocity profiles, 
a momentum thickness, 6 , of 2.55mm was obtained at x/H « —6 so that, for 
the flow at this position, Ro = 1680. 

The results from the three methods were: 

a) Cj = 3.854 x 10" 3 

b) Cj = 3.925 x 10" 3 

c) C f = 3.810 x 10" 3 

A mean value of Cj = 3.86 x 10 -3 was therefore used. This gave u T = 0.44m/s 
for a free stream velocity of lO.Om/s. The resulting velocity profiles, normalized 
using wall values, are shown in Figure 9. 

The curves are in good agreement except for the DISA probe. All curves have 
not been shown in Figures 8 and 9, but it was clear that the curves obtained 
with and without step were practically identical. The best fit straight line in 
the log region (Figure 9, neglecting the DISA results) is: 

U + = ~?—lny + + 5.35 
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1/2 

FIGURE 10. Distribution of u 2 and v 2 across the boundary layer at 

x/H « —6 (Outer scaling). 

1/2 

u 2 : • Miniature X-probe with step; o Miniature X-probe without step; 

A Single hot wire; x DISA X-probe 
1 /2 

v 2 : □ Miniature X-probe with step; V Miniature X-probe without step; 

• DISA X-probe. 

which is close to the fully developed case. 

The strength of the wake, using the miniature X-probe results, is A£/ + = 1.4. 
For a zero pressure gradient, Cole’s (1962) curve of n vs Re gives a wake strength 
of 1.90 for Re = 1680. Our wake strength indicates a slightly favorable pressure 
gradient. 

Based on the mean velocity profiles, it would appear that at 6 step heights 
upstream of the step, the flow is close to the fully developed state and that there 
is no influence of the step on the flow. 

9.4. Turbulence Quantities in Boundary Layer Upstream of Step 

Turbulence quantities from the single hot wire, from the miniature X-probe 
and from the DISA X-probe, obtained at the six step heights upstream of the 
step position, were compared with one another with published data. 

The u 2 results are shown in Figure 10 using outer scaling and in Figure 1 1 
using inner scaling. 
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x/H « —6 (Inner scaling). Symbols as for Figure 10. 

l fa 

u ? : — Spalart (1988); ••• Eckelmann (1974); x x x Andreopoulos et al. 

(1984). 

1/2 

v 2 : Spalart (1988); Eckelmann (1974); Andreopoulos et 

al. (1984) 


1/2 

All values in Figures 10 and 11 appear to be reasonable except for the v 2 
values obtained with the miniature X-probe. For this probe, the values obtained 
between y + = 10 to Y + = 100 (y = .36 to 3.6mm) are much higher than the 
usual values reported in the literature. Perry et al (1987) reported on and dis- 
cussed at some length the difficulty of obtaining consistent v 2 measurements in 
the inner regions of boundary layers. Their survey of 13 published results for 
smooth-wall boundary layers shows the scatter of results obtained and is repro- 
duced in Figure 12, together with the value that we obtained at that position in 
the boundary layer. 

Our maximum value agrees with the values reported at higher von Kasman 
numbers. These discrepancies need further investigation for clarification, al- 
though it should be _noted that the “non-isotropy,” indicated by the difference 
between the u 2 and v 2 values, seems rather high for the curves drawn in Figure 
11 . 

The uv results are shown in Figure 13 using outer scaling and in Figure 14 
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FIGURE 12. Values of v 2 reported in the literature at y/6 = 0.1 in a smooth- 
wall turbulent boundary layer. [Perry et al (1987)]. 
o Literature values — see Perry et al (1987) for details; 

• Value obtained with miniature X-probe at Y/6 = .1. 
x Maximum value obtained with miniature X-probe. 



FIGURE 13. Distribution of — uv across the boundary layer at x/H « —6 
(Outer scaling). 

• Miniature X Probe with step; o Miniature X-probe without step; x DISA 
X-probe. 


using inner scaling. 
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FIGURE 14. Distribution of — ut7 across the boundary layer at x/H « —6 
(Inner scaling). 

• Miniature X-probe with stop; o Miniature X-probe without step; x DISA 
X-probe. 

— Spalart (1988); • • • Eckelmann (1974); Andreopoulos et al (1984). 

Except for the Andreopoulos et al (1984) curve in the outer regions of the 
flow, all results show reasonable agreement for the utJ values. It is perhaps a 

little strange that results which are in agreement in u 2 (Figure 11) and uv 
(Figure 14) should be so different in the v 2 values (Figure 11). 

9.5. Characteristics of Tunnel Flow Approaching Close to Step 

When the Preston tube was traversed in the X direction along the plate (at 
z = 252mm — our standard z position) from x = —350mm to the lip of the 
step, we were surprised to find significant departures of Cj from the expected 
values — see Figure 15. 

This seemed to indicate an upstream influence of the step and, from pub- 
lished work on backward facing step flows, such an influence had not been 
expected. Consequently, using the Cj curve as a guide, five stations {x/h = 
—2.6, —1.5, —0.53, —0.13, and —0.) were selected for detailed profile studies. 
Since the miniature X-probe had given reasonable results, it was used for these 
profiles in order to allow measurements to be made as close to the wall as pos- 
sible. 
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arrows indicate stations where detailed profiles were taken. 



FIGURE 16. Distribution of uv in the outer regions of the boundary layer, 
upstream from the step. 

The changes indicated by the Cj profile (Figure 15) were reflected in the 
turbulence quantities - in particular in the uv values, as shown in Figures 16 
and 17. 

The minimum of —uv in the inner part of the flow, Figure 17, occurs at 
approximately the same point (x/H = —1.4) as the minimum in the C/ curve, 
Figure 15. In the outer part of the flow, Figure 16, there are obviously some 
complex changes occurring in the large scale structures. The negative values of 
—uv that were measured at y + = 10 cannot be substantiated at this stage and 
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FIGURE 17. Distribution of uv in the inner regions of the boundary layer, 
upstream from the step. 


may be due to the problems associated with X-wire measurements so close to 
the wall. A probe with wires of say 0.63/xm diameter and 0.1mm long may be 
required to resolve this region. 

The mean v velocities also follow the changes in the Cf curve, although in the 
opposite sense - peaking at minimum Cf and reaching a minimum at maximum 
Cj, Figure 18. 

This upstream influence of the step does not seem to have been previously 
reported in the literature, although three instances of support for our results 
have been obtained: 

a) In discussions with Dave Driver, he said that he had observed, in his work 
on the backward facing step, an increase in Cf (one data point) as the step 
was approached. This point was shown in the paper of Monson et al (1981) 
but was not commented on. The Cf values for their flow were different to 
ours, but by shifting the Cf origin and using the boundary layer thickness for 
non-dimensionalizing the distance from the step, some agreement is obtained 
between their Cf values and ours - see Figure 19. 

It is interesting to note that Monson et al (1981) obtained their Cf values 
using laser interferometry - a completely different technique to the Preston 
tube approach that we used. 

b) In 1981, M. Sindir carried out, for NASA-Ames, numerical simulation work, 
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FIGURE 18. Distribution of U in the inner regions of the boundary layer, 
upstream from the step. 



FIGURE 19. Cj values approaching the step. 

— Our results, Figure 15. x Monsom et al (1981) results. 

using modified TEACH programs, on flow over a backward facing step. In 
discussions with John Viegas concerning this work, John showed me the out- 
put from a program that he had run, using Sindir’s algebraic stress model 
program. It was observed that the pressure values, and these values only, did 
indicate some influence upstream of the step - the negative pressures in the 
recirculating zone extending upstream about two step heights from the step. 
The extent of the negative zone, above the surface, was about half a step 


I 


! 
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height at the step, decreasing to zero at 2.5 step heights from the step. These 
pressure values, however, show non-monotonic changes, indicating that a finer 
grid is required in this region of the flow, 
c) In his computational solution of the laminar flow over an aerofoil section, Jim 
Brown (EFD branch, NASA-Ames) showed that the Cj distribution had a 
distinct dip and then recovery just before reaching the trailing edge, in many 
ways emulating the curve of Figure 15. 

We conclude that the upstream influence of the step is a feature of the flow 
and may be an important consideration when direct simulations of the flow are 
attempted. 

9.6. Measurements Made Relating to the Large Scale Structures 

To study the large scale structures in the flow, it was decided to carry out 
initial measurements in the boundary layer upstream of the step. There were 
two reasons for this choice. Firstly, there is available considerable published 
material on the large scale structures in a boundary layer. Our results could 
then be of a corroborative nature and also perhaps, by judicious placement of 
the sensors and analysis of the data, contribute to the understanding of these 
structures - particularly with regard to the influence, one on the other, of the 
outer and inner structures. 

Secondly, any experimental and analysis techniques developed for this part of 
the flow could be applied equally well to the other regions of the flow — the 
“mixing layer” region and the region downstream of reattachment. All initial 
experiments were, therefore, carried out at x = —100mm (x/H = —2.63) where 
the step influence has hardly yet begun to have any major effect. For those 
measurements the floor plate was heated to a uniform temperature of about 
10° C above the free-stream temperature, as described earlier. 

The first experimental set-up consisted of two cold wire sensors at y + « 
10 from the wall and separated from each other by A z + « 100. The 3-wire 
miniature probe (X-wire + cold wire) was traversed between these two wires 
from y + « 10 to the outer regions of the boundary layer — data being taken at 
30 stations. A typical time trace from the five sensors is shown in Figure 20. 

Some structures can be discerned from the traces of the two wires closest to 
the wall. Sudden changes in the u signal also correspond to sudden changes in 
the temperature signals indicating that, as others have found, the u signal near 
the wall is a suitable structure detector. 

The second experimental set-up, for investigating the large scale structures, 
consisted of 14 sensors. The previous 5 sensors were complemented by a further 
5 cold wires arranged in a vertical array and 4 cold wires arranged in a horizontal 
array. Each of the wires in the vertical array was separated by approximately 
5mm (Ay + « 140) while the wire that was closest to the wall was at y + sa 90. 
To allow for the y traverse of the 3-wire sensor, it was necessary to offset the 
vertical array by about 4mm (A z + « 110) from the centre-line of the traverse. 
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FIGURE 20. Time traces of velocity and temperature fluctuations. 
x = —100mm, y + = 14. 0: on 3-wire probe, 0i and 02 : at y + = 10 


The horizontal array had 2 wires on one side of the traverse centre-line and 2 
wires on the other side. This array was situated at about the half-boundary-layer 
distance from the surface and spanned a total width of about 30mm (A z + « 
840). Data was obtained from all sensors as the 3-wire probe was moved to 30 
stations across the boundary layer. Typical traces from the 14 sensors, obtained 
when the 3-wire probe was at y = 14.8mm (y + = 415) are shown in Figure 21. 

A large scale structure can be noted as being detected simultaneously on a 
number of wires, as indicated by the crosses in Figure 21. The wire for 6j is 
obviously in the intermittent region of the flow. 

A program was written to enable large scale structures to be detected, using 
data similar to Figure 20 or 21, by any single channel or combination of channels, 
using the WAG, Bisset et al (1958), method. This program was run and the 
detections made compared favorably with the “by-eye” detections using the time 
traces. 

There was no time for further work by the author at this stage. 
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FIGURE 2 1 . Typical time traces of velocity and temperature fluctuations using 
3-wire probe ( u,v,0 ), 2 near-wall cold wires (0i, 02), a rake of 4 horizontal cold 
wires (03 to $q), and a rake of 5 vertical cold wires ( Oj to 0^). X = —100mm, 
y + = 415. 

10. Next Phase of the Work 

Broadly, the next phase of the work is to apply the detections made using the 
above program to obtain conditionally sampled values of the velocities (u and v) 
at each station across the boundary layer. These values should be obtained for 
the inner and outer regions of the flow and also for, say, only those structures 
with the most common spacing. The results for structures with other spacings 
will also need to be checked. Some correlation work is necessary to determine the 
relative position of the conditioned velocities. Topologies of the inner and outer 
regions can then be built up and the distribution (“coherent” and “random”) of 
momentum and heat transfer quantities determined. The other regions of the 
flow — the mixing layer and boundary layer after reattachment — can then be 
examined in a similar manner. 
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The effects of stabilizing and destabilizing 
longitudinal curvature on the structure 
of turbulent, two-stream mixing layers 

By M. W. PLESNIAK AND J. P. JOHNSTON 
1. Introduction 
1.1 Background 

The topic of the research program is mixing and turbulence transport in par- 
allel, co-flowing, two-stream shear or mixing layers. In particular the effects 
of streamwise (longitudinal) curvature on the growth and development of the 
shear layer were investigated. The preliminary results of the present research 
program are reported in Plesniak and Johnston (1987), where results are pre- 
sented for stable and unstable (in the Taylor-Gortler (T-G) sense) mildly curved 
two-stream mixing layers. These results are compared to a plane mixing layer 
studied in the same facility under nearly identical initial and boundary condi- 
tions. At typical operating conditions, a nominal speed ratio of 0.5 with flow 
speeds of 12 and 24 cm/s was employed. The ratio of the shear layer thickness 
to the radius of curvature is less than 5%. Reynolds numbers based on AU and 
momentum thickness ranged from 1000 to 2000 (5000 to 10000 based on vortic- 
ity thickness) over the region in which measurements were made. Laminar flow 
exists at the splitter plate trailing edge due to the low operating velocities. The 
average momentum thickness Reynolds number in the trailing edge boundary 
layers is estimated to be 190, too low for a trip to be effective. The mixing 
layer was not artificially forced and often exhibited more than one frequency 
of roll up. Dye flow visualization along with single component LDV results are 
presented for all three cases. 

Flow visualization indicated that for stabilizing curvature vortex mergings 
(pairings and triple interactions) are common, but there is little mixing of the 
dye across the layer and fine-scale three-dimensionality does not occur in the 
region studied. Conversely, destabilizing curvature promotes mixing at all scales 
and the layer is characterized by complex three-dimensional structures. Vortex 
interactions are also less prevalent for the unstable case than for the stable or 
plane layers. 

The LDV measurements showed that the unstable layer exhibits a greater 
growth than the plane or the stable layer of nearly the same velocity ratio 
(growth rate was determined from the maximum slope thickness). The mean 
(U) and fluctuating (u') streamwise velocity profiles exhibited self-similarity for 
all three cases studied. However the fluctuating normal («') profiles were not 
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self-similar for either the stable or unstable layer. Peak fluctuation levels in both 
components were greatest for the unstable layer. 

1.2 Current Research 

The main goal of the research program is to gain a fundamental understanding 
of the physics of curved, turbulent, unforced free mixing layers and to provide a 
data base for turbulence modelers. Since the facility was previously instrumented 
only with a one-component LDV system, the major objective was to construct a 
two-component LDV system to allow the measurement of Reynolds stresses and 
examine turbulence transport in the mixing layers. 

2. Approach 

2.1 Documentation of Initial Conditions 

It is well known that the scatter in data from various sources regarding free 
shear layers is large. This is due to the mixing layer’s extreme sensitivity to 
initial conditions. A number of studies have been undertaken to examine the 
effects of initial conditions and side-wall confinement effects on the growth of 
free shear layers. These studies attempted to resolve the discrepancies in the 
literature concerning spreading rates and other characteristics of mixing layers. 
However conflicts regarding the effects of laminar vs. turbulent boundary layers 
at separation, initial fluctuation level, initial momentum thickness and influence 
of splitter-plate trailing-edge thickness on the streamwise evolution of the shear 
layer still exist. Some groups report confinement effects of boundaries while 
others do not. 

In view of many unresolved issues concerning initial and boundary condi- 
tion effects it is considered important to document the initial conditions of the 
present study. The two-component LDV system described below is being used 
to investigate the boundary layers on the high and low speed sides of the splitter 
plate. Determination of the momentum thicknesses will allow us to predict the 
most unstable frequency and compare to the results in Ho and Huerre (1984). 
The levels of free stream turbulence in the nozzle will also be measured with 
this system. 

Another issue which has recently been addressed in the literature is spanwise 
nonuniformity in the mixing layer and the formation of small-scale streamwise 
vortical structures due to natural disturbances. Jimenez (1983), Lasheras et al 
(1986, 1988), Bernal and Roshko (1986), Mehta; Ho (private communications), 
and others have reported on the presence of these structures in plane mixing 
layers. They are easily detectable in spanwise variations of the mean streamwise 
(U) velocity profiles. Since these structures have never been studied in curved 
mixing layers, it is unknown whether there is any interaction with the T-G 
structures resulting in their suppression or enhancement. 
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2.2 Flow Visualization Studies 

Previous results obtained with vegetable dye visualization were able to show 
diffusion of the dye across the layer and also the entrainment of ambient fluid 
into the layer. VHS format videotape movies showed large-scale structures which 
transported dye across to opposite sides of the shear layer in the case of desta- 
bilizing curvature. However the dye visualization technique had several limita- 
tions: the major being the rapid diffusion of dye in the fully turbulent regimes 
restricting the operating velocities lower than those used in acquiring the quan- 
titative data. 

Thus, for the present study a hydrogen bubble technique with portable laser 
sheet illumination was employed. The wire may be positioned at any location 
in the mixing layer along the entire length of the test section. Using existing 
equipment, timelines can be generated by pulsing the wire. This technique will 
enable us to visualize the transport of turbulence (and of mass if the Schmidt 
number of the bubbles is low enough) by the structures of various scales. In 
particular, the large scale motions associated with the T-G rollers in the unstable 
layer can be studied. A suitable lighting source is currently being sought to 
illuminate a 2-D grid of hydrogen bubble wire used to visualize an entire plane 
across the mixing layer. 

2.3 Velocity Data 

After unexpected delays in delivery of the necessary components the two/three- 
component LDV system has been constructed, optically aligned and calibrated. 
Preliminary two-component velocity profiles have been obtained for the sta- 
ble and unstable layers. However, these data suffer convergence problems in the 
mixed (ur) moments due to insufficient numbers of samples. More complete data 
being presently acquired correct this problem. In addition the third component 
capability (w measured with an immersible fiber optics probe) has recently been 
added to the system so that the entire 3-D velocity field may simultaneously be 
acquired and the full Reynolds stress tensor computed. 

2.3.1 Description of the 2-Component LDV 

A three-color, two-component LDV system powered by a 3- Watt Argon-ion 
Lexel laser was constructed for turbulence studies in the curved mixing layers. 
The system employs frequency shifting on both channels and is configured to 
operate in the backscatter mode. The entire system is mounted on an opti- 
cal breadboard which may be traversed in a rectilinear X-Y global coordinate 
system. 

An aluminum traversing cart made of 2 x 4 inch square tubing (l/8 inch wall 
thickness) chosen for its light weight and stiffness was constructed to traverse 
the optical system. The traverse rides on pillow block linear bearings over a 1 
inch diameter ground stainless steel shaft. Positioning in the global streamwise 
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direction (X) is performed manually over a travel of 12 feet (3.66 m). The 
direction normal to this (V) has a travel of 4 feet (1.22 m). Thus, any point 
across the entire outer channel may be accessed. 

In order to traverse along a radial line in the inner curved test section, pro- 
grammable computer controlled stepper motors are used. The stepper motors 
are driven by Superior 430-TH translator modules operating in a chopped mode. 
An Anaheim Automation CL1605 programmable controller board is used to con- 
trol the translator modules. Parameters such as motor speed, acceleration, and 
limit switch condition are controlled by the board. A serial RS-232 interface is 
used to communicate to the CL1605 with the data acquisition computer system. 

A ball nut/screw assembly with a lead of 0.250 in (6.35 mm) is used for 
positioning in the Y direction while a lead screw unislide assembly with a lead 
of 0.100 in (2.54 mm) is employed for fine positioning over 10 inches (0.254 m) 
of travel in the global X direction. The stepper motors are driven in a half- 
step mode which produces a resolution of 400 steps/revolution. This allows a 
positioning resolution of 0.00025 in/step (0.00635 mm/step) in the X direction 
and 0.000625 in/step (0.015875 mm/step) in the Y direction. Backlash in X 
was measured to be 0.001 in (0.0254 mm) and in Y, 0.003 in (0.0762 mm). 

A radial move in the local x-y coordinate system of the curved shear layer 
consists of two moves in the global X-Y coordinate system. The angle between 
a radial line in the local coordinates and the global coordinate system is deter- 
mined experimentally. Considering spatial resolution in X and Y, backlash and 
uncertainty in determining the angle, the resolution in positioning accuracy is 
less than half of the LDV measuring volume diameter (100 microns). 

One limitation of this traversing system is that no automated spanwise travers- 
ing is possible due to design constraints resulting from lack of funding. However, 
the He-Ne based fiber optics probe traverse does have spanwise capability and 
can be used when such measurements are dictated. 

2.3.2 Present LDV Measurements 

The two-component LDV system described above is being used to make mea- 
surements in the curved mixing layer for three different configurations: straight 
(reference), stable and unstable with a velocity ratio of 2:1. Measurements of 
complete profiles (20 to 30 spatial points) at approximately 8 streamwise sta- 
tions will allow accurate determination of the mean growth rate and evolution 
of the shear layers’ turbulence fields. In conjunction with the fiber optics probe, 
three components of velocity (u,v,w) can be measured. Acquisition of all three 
velocity components simultaneously allows the entire Reynolds stress tensor to 
be computed and estimates of the turbulent energy budget can be made. Com- 
plete time records are saved at each position for post-processing. Mean and rms 
velocities, Reynolds stresses, higher order moments (skewness and kurtosis), and 
velocity probability density functions will be calculated. 
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Separate runs with much longer record lengths necessary to resolve the lowest 
frequencies will be made for spectral and autocorrelation data acquisition at 
several points in the mixing layer. Since we have two independent LDV systems 
(l-component probe and 2-component system), it will be possible to obtain two- 
point spatial correlations by positioning the probes at different points in the 
mixing layer. This kind of data may be particularly useful in tracking the large 
T-G rollers in the unstable layers. 

3. Summary 

A major portion of the period covered under CTR funding was spent on 
the construction and development of the multi-component traversing system 
and associated control hardware and software. Delays in delivery of necessary 
components delayed construction by several months. The system was first used 
to make measurements in October, 1988. 

In addition to acquiring the two-component measurement capability, a hydro- 
gen bubble/laser sheet flow visualization technique was developed to visually 
study the characteristics of the mixing layers. With this new technique we can 
look for evidence of the large-scale rollers arising from the Taylor-Gortler insta- 
bility and study its interaction with the primary Kelvin-Helmholtz structures. 

In general, the recently acquired data confirm the results reported in our 
preliminary paper (Plesniak and Johnston, 1987). At this stage we are making 
two and three-component measurements and will be able to report information 
regarding the transport of turbulence and full Reynolds stress data in several 
months. 
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An experimental investigation of a low 
Reynolds number turbulent boundary layer 
subject to an adverse pressure gradient 

By J. H. WATMUFF 


Summary 

This report covers a 12 month period from the beginning of the project. The 
project concerns an experimental study of a very low Reynolds number turbulent 
boundary layer subject to an adverse pressure gradient. The work is being per- 
formed in the Boundary Layer Wind Tunnel in the Fluid Mechanics Laboratory 
at NASA Ames Research Center in collaboration with Dr. R.V. Westphal. The 
aim is to obtain highly accurate mean-flow and turbulence measurements under 
conditions that can be closely related to the numerical simulations of Philippe 
Spalart for the purposes of CFD validation. It is expected that the experimental 
results will also serve as a useful contribution to the research literature in their 
own right. 

Much of the Boundary Layer Wind Tunnel has been completely rebuilt with 
a new wider contraction and working section which will improve compatibility 
with the simulations. A unique sophisticated high-speed computer controlled 3D 
probe traversing mechanism has also been integrated into the new test section. 
Construction of the tunnel and traverse has consumed a large fraction of the 
period covered by this report and is discussed in some detail. The hardware is 
now complete, and measurements are in progress. The mean-flow data indicate 
that a suitably two-dimensional base flow has been established. 

Automation of the probe positioning and data acquisition have led to a de- 
creased running time for total pressure measurements. However, the most sig- 
nificant benefits are expected to occur when using hot-wire probes. Calibrations 
can be performed automatically and there is no need to handle fragile probes 
when moving between measuring stations. Techniques are being developed which 
require sampling of the signals from moving hot-wire probes on the basis of their 
position in the flow. Measurements can be made in high intensity turbulence 
by flying probes upstream at high speed so that the relative magnitude of the 
turbulent velocity fluctuations are reduced. In regions where the turbulence in- 
tensity is not too large, the probe can also be repetitively scanned across very 
dense spatial grids in other directions. With this new technique, a complete 
profile can be measured in about 1/3 the time and with a spatial density about 
50 times that obtainable using a stationary probe. 
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1. Background 

1.1 CFD Simulations 

The numerical method used by Spalart (1988) is spectral in space and second- 
order in time. An important feature of the method is that there is no turbulence 
modelling. However, this requires a high grid density which restricts the simula- 
tions to low Reynolds numbers i.e. Re < 1500. Periodic spanwise and streamwise 
conditions are used and a multiple-scale procedure is applied to approximate slow 
streamwise growth. The key assumptions are that the streamwise evolution of 
the flow is slow and that the straining of the turbulence by the mean-flow can 
be neglected. Spalart has suggested that the assumption of small streamwise 
growth will cause the method to breakdown for large adverse pressure gradients. 

1.2 Relationship between the experiment and simulation 

There are three major requirements in the relationship between the experiment 
and simulation. Firstly, the Reynolds number of the experiment must match 
that of the simulation precisely. Secondly, the experiment and simulation also 
must have closely matched initial conditions. The simulation begins at the 
first station with an equilibrium boundary layer. A mildly favorable pressure 
gradient can be used to very closely approximate a “self-preserving layer” in 
the experiment, i.e. by careful experimental design the boundary layer can be 
maintained at almost constant thickness over some streamwise distance before 
being subjected to the adverse pressure gradient. Following a suggestion by 
Inman and Bradshaw (1981) the length of this region could be increased to allow 
any upstream trip effects to decay before the region of interest. Finally, accurate 
experimental pressure coefficient (C p ) measurements with high spatial resolution 
are required as an input for the simulation. A suitable flow configuration for 
the computations would be one in which the boundary layer experienced a non- 
dimensional pressure gradient ^ « 2 at a maximum Re « 1500. 

1.3 Preliminary experimental measurements 

The combined requirements of low Reynolds number and accurate pressure 
measurements precludes the use of water as a flow medium or the use of very low 
air velocities (e.g. 0.5 m/s) for the experiment. R.V. Westphal used momentum 
integral boundary layer calculations as a design tool to show that the Boundary 
Layer Wind Tunnel in its old form was a suitable facility for the experiment. 
Preliminary mean-flow measurements verified that this facility could duplicate 
the Reynolds number and pressure gradient requirements of the simulations. 
A flat test surface was employed for the boundary layer development, and the 
pressure gradient was imposed by a contoured upper control wall. The inlet flow 
velocity was 8.2 m/s. An inlet free-stream velocity around 8 m/s was considered 
to be close to the lower limit for accurate C p measurements, i.e. a total head of 
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0.15 inches of water with an accuracy of ±1%. The mean-flow two-dimensionality 
appeared to be good since the streamwise momentum balance was found to be 
within 10%. 

In the region of adverse pressure gradient at an Re « 1500 the reduced profile 
data showed a fl x « 0.7 which is a little low. A stronger adverse pressure gradient 
would be required for the current experiment. Also the inlet unit Reynolds 
number appeared to be a little high. A high inlet unit Reynolds number is 
undesirable since it shortens the physical distance over which the experiment 
matches the capabilities of the simulation. 

2. New Tunnel Contraction, Test Section and Traverse 

2.1 Tunnel Contraction 

When formulating the current experiment, it was desirable to use a lower inlet 
free-stream velocity which would reduce the inlet unit Reynolds number and also 
provide for better spatial and temporal resolution of hot-wire measurements. 
However, two factors had to be considered before a lower inlet velocity could 
be adopted. Firstly, the accuracy of the C p measurements would need to be 
maintained. Secondly, an increase in the total free-stream unsteadiness (from 
about 0.25% to around 0.5%) was observed at low tunnel speeds. 

The first restraint was relieved by purchasing the most accurate and sensitive 
commercially available pressure transducer. This allows the inlet free-stream 
velocity to be lowered to around 6 m/s (i.e. a total head around 0.08 inches of 
water) while still maintaining an accuracy of ±1% for the C p measurements. The 
second problem appeared to be related to the tunnel fan which was operating 
near the lower extremity of its characteristics. It was evident that a further drop 
in the fan speed using the original form of contraction and test section would 
lead to even more unsteadiness. Instead the flow area of the test section inlet 
was increased by 50% from 80cm x 20cm to 100cm x 24cm, allowing a lower 
inlet velocity but at a higher fan speed. This required designing a new 5:1 2D 
contraction to replace the original 7.5:1 3D version. 

Theoretical methods were used for selecting a new Boundary Layer Wind Tun- 
nel contraction design from a number of alternatives. The pressure distribution 
on the side walls was estimated using a potential flow method. The numerical 
solution of the stream function inside the contraction and inside upstream and 
downstream extensions of constant area was used to calculate the velocities on 
the boundary and the pressure distribution was obtained from the Bernouilli 
equation. The evolution of the contraction side wall laminar boundary layer 
subject to this pressure distribution was calculated using a momentum integral 
method. The effects of streamline curvature and convergence were neglected. 
The calculations showed that the new contraction design is not likely to cause 
separation. 
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2.2 High-speed probe traverse integral with test-section 

The numerical simulations usually have very high spatial resolution (e.g. 
Spalart (1988) used up to 10 7 grid points) but only for a limited number of 
time steps. This is in complete contrast to the experimental situation where 
almost unlimited temporal information is available but with very limited spatial 
content. One way of reducing this discrepancy is for the experimentalist to in- 
crease the spatial content of measurements. This dictates that multiple probes 
and/or more sophisticated traversing equipment be used. 

In view of the above considerations, a high-speed 3D computer-controlled 
probe traversing system has been built. The range of motion is 2.5m in the 
streamwise X-direction, 0.1m in the Y-direction normal to the wall and 0.5m 
in the spanwise Z-direction. Minimum top speeds of around 3 m/s in the X- 
direction and 1.5 m/s in the Y- and Z-directions were desired. These speeds 
can be reached over a short enough distance with accelerations of 2 g’s. Linear 
stepping motors satisfy the requirements for the Y- and Z-axes for accurate (0.05 
mm) high speed (1.5 m/s) positioning. The probe and sting are carried directly 
by the Y-axis motor. The Y-axis motor rides on a steel platen which is attached 
to the Z-axis motor, linear bearing and platen assembly. The Z-axis assembly 
is supported above the test wall within the working section by a gantry which 
spans the lm distance between the side walls. 

The requirements of the gantry are low mass (to minimize inertial loading), 
high stiffness (to minimize deflections), and small projected area (to minimize 
aerodynamic interference). Therefore, the gantry was constructed in the form 
of a lm wide truss using sections cut from thin sheets of carbon-fiber composite 
which were bonded together with epoxy resin. By itself the gantry weighs about 
2Kg which is less than 10% of the total movable mass. The total projected 
area of everything (i.e. the gantry, bearing rails, motors, platens etc.) located 
within the tunnel is equal to a one inch thick object spanning the tunnel side 
walls so the blockage is not excessive. Loading the gantry with a lOKg weight 
in the center of its one meter span results in a deflection of only 0.05mm. The 
ends of the gantry are fixed to carriages which move on linear bearings in the 
X-direction. 

The gantry is propelled in the X-direction by a brushless linear d.c. motor 
which reacts with a long magnet track to provide a maximum propulsive force 
of 80lbf. The motor is operated as a closed loop servo system with feedback 
provided by a linear quadrature encoder with a resolution of 10 micron. The 
maximum speed of the gantry has been tested at 3 m/s. Current linear stepping 
motor technology cannot match this performance. The two linear bearing rails 
are fixed to adapter plates which are bolted to a 1.2 x 3.0 m (4 x 10 foot) optics 
table. One of these adapter plates also supports the magnet track and linear 
encoder. The optics table provides an extremely rigid mounting platform for the 
entire assembly. Another significant benefit of the optics table is that it provides 
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an extremely flat reference surface for alignment purposes. 

Measurements are made in the boundary layer forming on a test plate which 
consists of a 0.375 ± 0.005 inch thick ground aluminum toolplate approximately 
lm wide by 2.1m in length. The plate is supported above the optics table by 
machined standoffs. A dial gauge attached to the traverse indicated that the 
deviations of the test-plate over the entire traversable X-Z plane were within 
±0.15 mm. The X-axis motor and twin-rail linear bearings are located beneath 
the test-plate, and rubber strips are attached to the side walls and to the plate 
to seal the slots that provide access for mounting the gantry. The side walls and 
flexible ceiling which complete the test-section are attached to a frame that is 
bolted to the optics table. Full access to the test-plate and traverse assembly 
can be obtained in minutes by releasing the bolts and hoisting the frame away 
with a crane. 

3. Benefits of the Traverse-Some New Techniques 

Automation of the high-speed 3D traverse reduces the time taken for experi- 
ments since repetitive control commands do not have to be performed manually. 
Although this leads to decreased running time for total pressure measurements, 
the most significant benefits occur when using hot-wire probes. 

Hot-wire probes can only be used with an acceptable accuracy when the mean 
velocity component is large compared to the velocity fluctuations. For a normal 
wire, the assumption of sensitivity only to the streamwise velocity fluctuations 
will become poorer as the fluctuations of the o'ther component normal to the 
wire increase (by simple vector addition). Crossed-wire probes have a theoret- 
ical velocity vector wedge-angle limit defined by the normals to each wire. In 
practice, the allowable flow angle is much less than this owing to wakes which are 
shed from the prongs. The boundary layer turbulence encountered in the region 
of adverse pressure gradient is expected to increase to such an extent that these 
sources of error will become significant. One way around this problem is to fly 
the probe in the upstream direction at high speed while sampling on the basis of 
the probe position. The superimposed bias velocity will reduce the magnitude 
of the turbulent velocity fluctuations relative to the increased velocity seen by 
the probe so that accurate measurements are possible. 

There is another, very distinct, advantage of flying the wire, i.e. measurements 
can be obtained on a very dense spatial grid. The time required to obtain 
converged data by repetitively scanning a dense grid can be more than an order 
of magnitude less than the time required to get the same data by traversing the 
probe to each point on the grid in succession and taking data with a stationary 
probe. At the low speeds of this experiment the sampling period must be as 
long as 100 seconds for adequate data convergence. In this case, the total time 
to obtain 30 points in a profile in the usual way with a stationary probe is 
of about 50 minutes duration. In regions where the turbulence intensity is 
not too high, the probe can be scanned up and down through the layer at 
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speeds that are quite high but not large enough relative to the mean flow to 
exceed the thresholds mentioned above. For example, with constant magnitude 
accelerations of alternating sign and a maximum speed limit of around 0.6 m/s, 
the oscillation frequency for traversing a distance of 32mm up and back down 
through a layer is about 4.5Hz. If samples are taken on the basis of the probe 
position and if 5,000 samples at each point are adequate for data convergence, 
then only 18 minutes is required to measure the entire profile compared to 50 
minutes in the case of the stationary probe. However, a more significant benefit 
is that there is no time penalty for increasing the spatial density of the grid. For 
example, if samples are taken for each 0.025mm increment of the probe position, 
then the spatial density of the grid can be increased by a factor of more than 40 
compared to the conventional method with a stationary probe. Also, spanwise 
derivatives and streamwise derivatives can be measured. 

Calibration of crossed-wire probes requires a mechanism in addition to the 
probe traverse. For example, in dynamic calibration schemes the small per- 
turbation sensitivity of the probe is determined directly by oscillating it back 
and forth continuously in a steady and uniform free-stream using a mechanical 
shaker. An alternative calibration scheme consists of imposing a number of flow 
angles on the probe by rotating it about an axis normal to the plane of the 
prongs. Intricate mechanisms are usually added to the traverse to perform this 
function. However, the same effect can be accomplished relative to the probe 
by moving it at high speed across a uniform free-stream. The advantage of a 
high speed 3D traverse is that either crossed-wire calibration scheme can be 
used without the need for additional hardware and without the need to handle 
probes. 

Hot-wire probes must be calibrated frequently in a uniform stream. Using a 
conventional single-axis traversing system in the adverse pressure gradient ex- 
periments would require the removal and reinstallation of fragile and expensive 
probes from the measurement region for the purpose of calibration. This would 
be time consuming and risk probe breakage. With a high-speed traversing sys- 
tem this could be performed almost instantly as well as minimizing the chances 
of probe damage. 

Since all the functions described above are under software control, the oppor- 
tunity exists to create an intelligent system. For example, hot-wire calibration 
drift is a particularly common and frustrating phenomenon during the course 
of an experiment. At regular intervals, the probe could be rapidly moved to 
a reference point in the free-stream to check for drift and another calibration 
performed (automatically) if some tolerance level was exceeded. 

4. Results to Date 

4.1 Incoming laminar boundary layer 

The boundary layer on the test plate has its origin upstream and, therefore, 
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it is influenced by the pressure gradients, wall curvature and lateral conver- 
gence within the contraction. Since the flow is accelerated by a contraction, 
the pressure gradients on the side walls are generally favourable. However, it 
can be shown that for a finite length contraction two regions of adverse pres- 
sure gradient exist on the walls in the vicinity of the entrance and the exit. 
Laminar boundary layers are especially sensitive to adverse pressure gradients 
and are prone to separation causing increased unsteadiness. Because of the low 
Reynolds number requirements of the current experiment, the quality of the 
incoming laminar boundary layer needs even more attention than usual. 

Laminar boundary layer profiles were measured with a total pressure probe 
at a station 8.3cm downstream of the contraction exit for a number of free- 
stream velocities. The pressure was constant along the test section for these 
measurements. All the velocity profiles were found to be very near the Blasius 
profile. In particular the boundary layer corresponding to the lowest velocity of 
U 0 = 6.3 m/s planned for the experiment has a shape factor H=^-=2.45 which 
is close to that of the Blasius profile (H=2.59). 

4.2 Selection of a transition device 

The low Reynolds number requirements of the experiment dictates that mea- 
surements be obtained close to the boundary layer tripping device. The effect 
of a simple cylindrical wire has been the subject of many studies, see Schlicht- 
ing (1979) p. 537 for a review. More recently there have been observations of 
spanwise irregularities in the boundary layers behind trip wires. Many workers 
have suggested that distributed three-dimensional roughness elements may be 
superior for transition purposes. However, the author is not aware of a system- 
atic parametric study that offers a reproducible alternative to a trip wire. The 
observed irregularities have varying strengths depending on the facility. There is 
evidence to suggest (see Bradshaw 1965) that the observations are more closely 
associated with wind tunnel screens rather than with trip wires. 

For the reasons outlined above, a simple wire was selected for boundary layer 
transition. The diameter and streamwise position for a number of wires were 
determined using the guidelines reported in Schlichting in conjunction with the 
laminar velocity profile measurements discussed earlier. Satisfactory results were 
obtained with 1.7mm and 2.0mm (diameter) wires when positioned near the 
contraction exit. However, when moved to a location 20cm downstream of the 
contraction, the 1.7mm wire failed to produce transition at all. Transition was 
intermittent with the 2.0mm wire at this location. These observations can be 
explained in terms of the adverse pressure gradient (not measured) near the 
exit of the contraction. Effective transition was obtained with a larger 2.4mm 
diameter wire at X=15cm. Measurements of Cj, momentum thickness and 
shape factor H indicated that a regular turbulent boundary layer is established 
behind the 2.4mm wire by X « 25cm corresponding to a momentum thickness 
Re « 450. 
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4.3 Mean velocity profiles with the pressure gradient 

The variation of C p with streamwise distance is shown in figure 1(a). The 
streamwise spacing (5cm) of the measurements appears to be adequate since 
estimates of dC p /dX obtained with four different numerical schemes are much 
the same as shown in figure 1(b). Mean velocity profiles were measured on the 
centerline (Z=0) with a 0.042 inch outside diameter total pressure probe which 
also functioned as a Preston tube for estimating the wall shear stress. The same 
probe was also used for the spanwise wall shear stress surveys. Mean pressures 
were determined over a 90 second period requiring about 45 minutes per profile 
and the data was acquired over a three day period. The probe was traversed 
under-computer control and once set running in the morning, the experiment 
proceeded without manual intervention throughout the day. 

A number of interesting features can be determined from the data set. A list 
of quantities derived from the profile data is given in table 1 and some of these 
have been plotted in figures 2(a)-(f). For 35cm < X < 55cm the layer thicknesses 
are approximately constant and the shape of the velocity profiles is similar as 
indicated by the shape factor. This suggests that the layer in this region is close 
to equilibrium. The profiles are shown in both wall- and outer-coordinates in 
both the region of the favorable and the adverse pressure gradients in figures 
3(a)-(d). Estimates of the skin friction coefficient Cf 0 obtained from the Preston 
tube have been used to plot the profiles shown in wall-coordinates. The pressure 
gradients are well within the allowable limits suggested by Patel (1965). The 
deviations from the log law indicate that uncertainties of up to 10% in wall shear 
stress may have to be accepted. 

4.4 Checks for two-dimensionality 

The results of two spanwise C/ surveys are shown in figure 4(a). The location 
of the upstream survey corresponds to the region where the layer is in equilib- 
rium. The survey spans approximately 40 boundary layer thicknesses and the 
variations are within ± 5% of the mean value. The downstream location in the 
region of adverse pressure gradient corresponds to Re « 1400 which is within the 
capabilities of Spalart’s simulations. The spanwise extent of the survey is about 
20 layer thicknesses, and again the variations are within ± 5% of the mean value. 
The momentum balance shown in figure 4(b) is within 10%. The spanwise Cj 
measurements and the momentum balance indicate that the flow is acceptably 
two-dimensional. 

5. Conclusions 

The boundary layer in the remodelled Boundary Layer Wind Tunnel is closely 
related to the capabilities of Spalart’s simulations. In particular, the layer ap- 
pears to be close to equilibrium just before application of the adverse pressure 
gradient which is a suitable initial condition for a simulation. The pressure gra- 
dient parameter /3 X « 2 at R e « 1650 which is about the upper limit of the 
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numerical simulation capabilities. At this location, the distance downstream of 
the contraction is equal to the width of the test section. The spanwise variation 
of Cf is within ± 5% and the momentum balance is good. 

The high-speed computer-controlled 3D traverse has unique capabilities that 
will allow hot-wire measurements with a high spatial density to be obtained very 
quickly. Accurate measurements can be obtained on the “fly” in regions of high 
turbulence intensity that would be otherwise be subject to gross uncertainties 
using stationary probes. 
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TAB LEI. Streamwise Development of Mean Velocity Profiles with Pressure Gradient 


NOTE: 

1. Reference inlet velocity U # %6.$ m/s at X=Qcm. 

2. Boundary layer tr : p wire k at X=15cxn. 
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Figure 1(a) Variation of Pressure Coefficient C ? with streamwise distance, 
(b) Estimates of dC p /dX using four numerical schemes. 
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Figure 2. Variation of quantites derived from profile data with streamwise distance. 

(a) Skin Friction Coefficient C/ (b) Displacement 6* and momentum 0 thicknesses 
(c) Reynolds number based on momentum thickness R $ (d) Boundary layer thickness 6 
(e) Shape Factor H=y. (f) Pressure gradient parameter 
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Figure 3. Mean velocity profiles. 

(a) Region of favorable pressure gradient and (L) in wall coordinates, 
(c) Region of adverse pressure gradient and (d) in wall coordinates. 





Figure 4(a) Spanwise variation of C / at X=50cm and X=90cm 
(b) Momentum balance. 
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APPENDIX 


CENTER FOR TURBULENCE RESEARCH 
1988 ROSTER 


NAME/TERM 


AREA OF RESEARCH 

POSTDOCTORAL FELLOWS 



CABOT, Dr. William H. 
3/88 - present 

Ph.D. Physics, 1983, University 
of Rochester 

Turbulence in the early 
solar nebula, convection 

KARNIADAKIS, Dr. George 
9/87 - 2/88 

Ph.D. Mechanical Engineering, 
Massachusetts Institute of 
Technology 

Simulation and analysis of 
basic complex geometry 
turbulent flows 

LEE, Dr. Moon J. 
12/87 - present 

Ph.D. Mechanical Engineering, 
1985, Stanford 

Turbulence physics and 
modeling 

LELE, Sanjiva 
1/88 - present 

Ph.D. Mechanical Engineering, 
1985, Cornell University 

Direct numerical 
simulation of 
compressible free shear 
flows 

POINSOT, Dr. Thierry 
9/88 - present 

Docteur es Sciences, 
Mechanical Engr., 1987, Univ. 
d’Orsay, France 

Direct simulation of 
turbulent reacting flows 

SHIH, Dr. Tsan Hsing 
4/87 - present 

Ph.D Mechanical Engineering, 
1984, Cornell University 

Turbulence modeling 

SMITH, Dr. Leslie M. 
9/88 - present 

Ph.D. Applied Mathematics, 
1988, Massachusetts Institute of 
Technology 

Renormalization group 
theory of turbulence 

STANAWAY, Dr. Sharon K. 
6/88 - 12/88 

Ph.D. Aeronautics and 
Astronautics, 1988, Stanford 

Numerical simulation of 
viscous vortex rings 

THOMPSON, Dr. Kevin W. 
6/87 - present 

Ph.D. Physics, 1985, Princeton 

Direct numerical 
simulation of 
compressible turbulent 
flows with application to 
the early evolution of solar 
nebula 

VASTANO, Dr. John A. 
9/88 - present 

Ph.D. Physics, 1988, University 
of Texas at Austin 

Low dimensional chaos in 
turbulence 
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NAME/TERM AREA OF RESEARCH 


VEGT, Dr. Jaap van der 
9/88 - present 

Ph.D. Mathematics, 1988, Delft 
Univ. of Technology, The 
Netherlands 

Transition to turbulence at 
hypersonic speeds 

WATMUFF, Dr. Jonathan H. 
1 1/87 - present 

Ph.D. Mechanical Engineering, 
1979, University of Melbourne, 
Australia 

Experimental 
investigation of turbulent 
boundary layers with 
adverse pressure gradient 

SENIOR POSTDOCTORAL FELLOWS 


BROADWELL, Dr. James E. 
8/88 - present 

California Institute of 
Technology 

Turbulent reacting flows 

BROWNE, Dr. L. W. B. 
2/88 - 8/88 

University of Newcastle - 
Australia 

Experimental 
investigation of flow over 
backward facing step 

ETEMADI, Prof. Nassrollah 
8/87 - 5/88 

University of Illinois 

Lagrangian statistics 

GERMANO, Prof. Massimo 
11/88 

University of Torino - Italy 

Subgrid scale modeling 
and filtering 

LESIEUR, Prof. Marcel 
4/88 - 5/88 

University of Grenoble, France 

Two-point closures 

LILLY, Prof. Douglas K. 
8/88 - 12/88 

CIMMS - University of 
Oklahoma 

Subgrid scale modeling 

MELANDER, Dr. Mogens 
5/88 - 7/88 

Univ. of Pittsburgh 

3D vortex dynamics at 
high Reynolds numbers 

ROSHKO, Prof. Anatol 
4/88 

California Institute of 
Technology 

Turbulent free shear flows 
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NAME/TERM 

AREA OF RESEARCH 

GRADUATE STUDENTS 

BEAUDAN, Patrick 
10/87 - present 

Numerical methods 

COLONIUS, Timothy 
7/88 - 9/88 

Aerodynamic noise 

KASSINOS, Stavros 

Reynolds averaged 

10/88 - present 

turbulence modeling 

LE, Hung 

Direct numerical 

4/88 - present 

simulation of flow over a 
backward facing step 

LEE, Sangsan 

Shock/turbulence 

4/88 - 12/88 

interaction 

NEVES, Joao 

Numerical simulation of 

10/87 - 4/88 

axial flow over a cylinder 

PLESNIAK, Michael 

Effects of curvature on the 

1/88 - 6/88 

structure of mixing layers 

POLIFKE, Wolfgang 

Topological complexity of 

City Univ. of NY 

the vorticity field in 

10/88 - 12/88 

homogeneous turbulence 

SASGES, Michael 

Mixing in curved shear 

7/88 - 9/88 

layers 

SORAKAYALPET, 

Space-time characteristics 

Sundaresan 

of wall -pressure 

10/87 - 4/88 

fluctuations 

SQUIRES, Kyle 

Effects of particle loading 

10/87 - 9/88 

on turbulent structure and 
modeling 



